1 /* mpn_sec_invert
2 
3    Contributed to the GNU project by Niels Möller
4 
5 Copyright 2013 Free Software Foundation, Inc.
6 
7 This file is part of the GNU MP Library.
8 
9 The GNU MP Library is free software; you can redistribute it and/or modify
10 it 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 
16 or
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 
22 or both in parallel, as here.
23 
24 The GNU MP Library is distributed in the hope that it will be useful, but
25 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
26 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
27 for more details.
28 
29 You should have received copies of the GNU General Public License and the
30 GNU Lesser General Public License along with the GNU MP Library.  If not,
31 see 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. */
38 static mp_size_t
39 mpn_cnd_neg_itch (mp_size_t n)
40 {
41   return n;
42 }
43 #endif
44 
45 /* FIXME: Ought to return carry */
46 static void
mpn_cnd_neg(int cnd,mp_limb_t * rp,const mp_limb_t * ap,mp_size_t n,mp_ptr scratch)47 mpn_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 
54 static int
mpn_sec_eq_ui(mp_srcptr ap,mp_size_t n,mp_limb_t b)55 mpn_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 
68 mp_size_t
mpn_sec_invert_itch(mp_size_t n)69 mpn_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. */
81 int
mpn_sec_invert(mp_ptr vp,mp_ptr ap,mp_srcptr mp,mp_size_t n,mp_bitcnt_t bit_size,mp_ptr scratch)82 mpn_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