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.h"
34 #include "gmp-impl.h"
35 
36 #if 0
37 /* Currently unused. Should be resurrected once mpn_cnd_neg is
38    advertised. */
39 static mp_size_t
40 mpn_cnd_neg_itch (mp_size_t n)
41 {
42   return n;
43 }
44 #endif
45 
46 /* FIXME: Ought to return carry */
47 static void
mpn_cnd_neg(int cnd,mp_limb_t * rp,const mp_limb_t * ap,mp_size_t n,mp_ptr scratch)48 mpn_cnd_neg (int cnd, mp_limb_t *rp, const mp_limb_t *ap, mp_size_t n,
49 	     mp_ptr scratch)
50 {
51   mpn_lshift (scratch, ap, n, 1);
52   mpn_cnd_sub_n (cnd, rp, ap, scratch, n);
53 }
54 
55 static void
mpn_cnd_swap(int cnd,volatile mp_limb_t * ap,volatile mp_limb_t * bp,mp_size_t n)56 mpn_cnd_swap (int cnd, volatile mp_limb_t *ap, volatile mp_limb_t *bp,
57 	      mp_size_t n)
58 {
59   volatile mp_limb_t mask = - (mp_limb_t) (cnd != 0);
60   mp_size_t i;
61   for (i = 0; i < n; i++)
62     {
63       mp_limb_t a, b, t;
64       a = ap[i];
65       b = bp[i];
66       t = (a ^ b) & mask;
67       ap[i] = a ^ t;
68       bp[i] = b ^ t;
69     }
70 }
71 
72 static int
mpn_sec_eq_ui(mp_srcptr ap,mp_size_t n,mp_limb_t b)73 mpn_sec_eq_ui (mp_srcptr ap, mp_size_t n, mp_limb_t b)
74 {
75   mp_limb_t d;
76   ASSERT (n > 0);
77 
78   d = ap[0] ^ b;
79 
80   while (--n > 0)
81     d |= ap[n];
82 
83   return d == 0;
84 }
85 
86 mp_size_t
mpn_sec_invert_itch(mp_size_t n)87 mpn_sec_invert_itch (mp_size_t n)
88 {
89   return 4*n;
90 }
91 
92 /* Compute V <-- A^{-1} (mod M), in data-independent time. M must be
93    odd. Returns 1 on success, and 0 on failure (i.e., if gcd (A, m) !=
94    1). Inputs and outputs of size n, and no overlap allowed. The {ap,
95    n} area is destroyed. For arbitrary inputs, bit_size should be
96    2*n*GMP_NUMB_BITS, but if A or M are known to be smaller, e.g., if
97    M = 2^521 - 1 and A < M, bit_size can be any bound on the sum of
98    the bit sizes of A and M. */
99 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)100 mpn_sec_invert (mp_ptr vp, mp_ptr ap, mp_srcptr mp,
101 		mp_size_t n, mp_bitcnt_t bit_size,
102 		mp_ptr scratch)
103 {
104   ASSERT (n > 0);
105   ASSERT (bit_size > 0);
106   ASSERT (mp[0] & 1);
107   ASSERT (! MPN_OVERLAP_P (ap, n, vp, n));
108 #define bp (scratch + n)
109 #define up (scratch + 2*n)
110 #define m1hp (scratch + 3*n)
111 
112   /* Maintain
113 
114        a = u * orig_a (mod m)
115        b = v * orig_a (mod m)
116 
117      and b odd at all times. Initially,
118 
119        a = a_orig, u = 1
120        b = m,      v = 0
121      */
122 
123 
124   up[0] = 1;
125   mpn_zero (up+1, n - 1);
126   mpn_copyi (bp, mp, n);
127   mpn_zero (vp, n);
128 
129   ASSERT_CARRY (mpn_rshift (m1hp, mp, n, 1));
130   ASSERT_NOCARRY (mpn_sec_add_1 (m1hp, m1hp, n, 1, scratch));
131 
132   while (bit_size-- > 0)
133     {
134       mp_limb_t odd, swap, cy;
135 
136       /* Always maintain b odd. The logic of the iteration is as
137 	 follows. For a, b:
138 
139 	   odd = a & 1
140 	   a -= odd * b
141 	   if (underflow from a-b)
142 	     {
143 	       b += a, assigns old a
144 	       a = B^n-a
145 	     }
146 
147 	   a /= 2
148 
149 	 For u, v:
150 
151 	   if (underflow from a - b)
152 	     swap u, v
153 	   u -= odd * v
154 	   if (underflow from u - v)
155 	     u += m
156 
157 	   u /= 2
158 	   if (a one bit was shifted out)
159 	     u += (m+1)/2
160 
161 	 As long as a > 0, the quantity
162 
163 	   (bitsize of a) + (bitsize of b)
164 
165 	 is reduced by at least one bit per iteration, hence after (bit_size of
166 	 orig_a) + (bit_size of m) - 1 iterations we surely have a = 0. Then b
167 	 = gcd(orig_a, m) and if b = 1 then also v = orig_a^{-1} (mod m).
168       */
169 
170       ASSERT (bp[0] & 1);
171       odd = ap[0] & 1;
172 
173       swap = mpn_cnd_sub_n (odd, ap, ap, bp, n);
174       mpn_cnd_add_n (swap, bp, bp, ap, n);
175       mpn_cnd_neg (swap, ap, ap, n, scratch);
176 
177       mpn_cnd_swap (swap, up, vp, n);
178       cy = mpn_cnd_sub_n (odd, up, up, vp, n);
179       cy -= mpn_cnd_add_n (cy, up, up, mp, n);
180       ASSERT (cy == 0);
181 
182       cy = mpn_rshift (ap, ap, n, 1);
183       ASSERT (cy == 0);
184       cy = mpn_rshift (up, up, n, 1);
185       cy = mpn_cnd_add_n (cy, up, up, m1hp, n);
186       ASSERT (cy == 0);
187     }
188   /* Should be all zeros, but check only extreme limbs */
189   ASSERT ( (ap[0] | ap[n-1]) == 0);
190   /* Check if indeed gcd == 1. */
191   return mpn_sec_eq_ui (bp, n, 1);
192 #undef bp
193 #undef up
194 #undef m1hp
195 }
196