1/* mpn_sec_invert
2
3   Contributed to the GNU project by Niels M��ller
4
5Copyright 2013 Free Software Foundation, Inc.
6
7This file is part of the GNU MP Library.
8
9The GNU MP Library is free software; you can redistribute it and/or modify
10it under the terms of either:
11
12  * the GNU Lesser General Public License as published by the Free
13    Software Foundation; either version 3 of the License, or (at your
14    option) any later version.
15
16or
17
18  * the GNU General Public License as published by the Free Software
19    Foundation; either version 2 of the License, or (at your option) any
20    later version.
21
22or both in parallel, as here.
23
24The GNU MP Library is distributed in the hope that it will be useful, but
25WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
26or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
27for more details.
28
29You should have received copies of the GNU General Public License and the
30GNU Lesser General Public License along with the GNU MP Library.  If not,
31see https://www.gnu.org/licenses/.  */
32
33#include "gmp-impl.h"
34
35#if 0
36/* Currently unused. Should be resurrected once mpn_cnd_neg is
37   advertised. */
38static mp_size_t
39mpn_cnd_neg_itch (mp_size_t n)
40{
41  return n;
42}
43#endif
44
45/* FIXME: Ought to return carry */
46static void
47mpn_cnd_neg (int cnd, mp_limb_t *rp, const mp_limb_t *ap, mp_size_t n,
48	     mp_ptr scratch)
49{
50  mpn_lshift (scratch, ap, n, 1);
51  mpn_cnd_sub_n (cnd, rp, ap, scratch, n);
52}
53
54static int
55mpn_sec_eq_ui (mp_srcptr ap, mp_size_t n, mp_limb_t b)
56{
57  mp_limb_t d;
58  ASSERT (n > 0);
59
60  d = ap[0] ^ b;
61
62  while (--n > 0)
63    d |= ap[n];
64
65  return d == 0;
66}
67
68mp_size_t
69mpn_sec_invert_itch (mp_size_t n)
70{
71  return 4*n;
72}
73
74/* Compute V <-- A^{-1} (mod M), in data-independent time. M must be
75   odd. Returns 1 on success, and 0 on failure (i.e., if gcd (A, m) !=
76   1). Inputs and outputs of size n, and no overlap allowed. The {ap,
77   n} area is destroyed. For arbitrary inputs, bit_size should be
78   2*n*GMP_NUMB_BITS, but if A or M are known to be smaller, e.g., if
79   M = 2^521 - 1 and A < M, bit_size can be any bound on the sum of
80   the bit sizes of A and M. */
81int
82mpn_sec_invert (mp_ptr vp, mp_ptr ap, mp_srcptr mp,
83		mp_size_t n, mp_bitcnt_t bit_size,
84		mp_ptr scratch)
85{
86  ASSERT (n > 0);
87  ASSERT (bit_size > 0);
88  ASSERT (mp[0] & 1);
89  ASSERT (! MPN_OVERLAP_P (ap, n, vp, n));
90#define bp (scratch + n)
91#define up (scratch + 2*n)
92#define m1hp (scratch + 3*n)
93
94  /* Maintain
95
96       a = u * orig_a (mod m)
97       b = v * orig_a (mod m)
98
99     and b odd at all times. Initially,
100
101       a = a_orig, u = 1
102       b = m,      v = 0
103     */
104
105
106  up[0] = 1;
107  mpn_zero (up+1, n - 1);
108  mpn_copyi (bp, mp, n);
109  mpn_zero (vp, n);
110
111  ASSERT_CARRY (mpn_rshift (m1hp, mp, n, 1));
112  ASSERT_NOCARRY (mpn_sec_add_1 (m1hp, m1hp, n, 1, scratch));
113
114  while (bit_size-- > 0)
115    {
116      mp_limb_t odd, swap, cy;
117
118      /* Always maintain b odd. The logic of the iteration is as
119	 follows. For a, b:
120
121	   odd = a & 1
122	   a -= odd * b
123	   if (underflow from a-b)
124	     {
125	       b += a, assigns old a
126	       a = B^n-a
127	     }
128
129	   a /= 2
130
131	 For u, v:
132
133	   if (underflow from a - b)
134	     swap u, v
135	   u -= odd * v
136	   if (underflow from u - v)
137	     u += m
138
139	   u /= 2
140	   if (a one bit was shifted out)
141	     u += (m+1)/2
142
143	 As long as a > 0, the quantity
144
145	   (bitsize of a) + (bitsize of b)
146
147	 is reduced by at least one bit per iteration, hence after (bit_size of
148	 orig_a) + (bit_size of m) - 1 iterations we surely have a = 0. Then b
149	 = gcd(orig_a, m) and if b = 1 then also v = orig_a^{-1} (mod m).
150      */
151
152      ASSERT (bp[0] & 1);
153      odd = ap[0] & 1;
154
155      swap = mpn_cnd_sub_n (odd, ap, ap, bp, n);
156      mpn_cnd_add_n (swap, bp, bp, ap, n);
157      mpn_cnd_neg (swap, ap, ap, n, scratch);
158
159      mpn_cnd_swap (swap, up, vp, n);
160      cy = mpn_cnd_sub_n (odd, up, up, vp, n);
161      cy -= mpn_cnd_add_n (cy, up, up, mp, n);
162      ASSERT (cy == 0);
163
164      cy = mpn_rshift (ap, ap, n, 1);
165      ASSERT (cy == 0);
166      cy = mpn_rshift (up, up, n, 1);
167      cy = mpn_cnd_add_n (cy, up, up, m1hp, n);
168      ASSERT (cy == 0);
169    }
170  /* Should be all zeros, but check only extreme limbs */
171  ASSERT ( (ap[0] | ap[n-1]) == 0);
172  /* Check if indeed gcd == 1. */
173  return mpn_sec_eq_ui (bp, n, 1);
174#undef bp
175#undef up
176#undef m1hp
177}
178