Line data Source code
1 : #include "fd_f25519.h"
2 : #include "../hex/fd_hex.h"
3 :
4 : #if FD_HAS_AVX512
5 : #include "avx512/fd_f25519.c"
6 : #else
7 : #include "ref/fd_f25519.c"
8 : #endif
9 :
10 : #ifdef FD_HAS_S2NBIGNUM
11 : #include <stdint.h>
12 : #include "../../third_party/s2n-bignum/include/s2n-bignum.h"
13 : #include "../../util/sanitize/fd_msan.h"
14 : #endif
15 :
16 : void
17 : fd_f25519_debug( char const * name,
18 0 : fd_f25519_t const * a ) {
19 0 : char *
20 0 : fd_hex_encode( char * FD_RESTRICT dst,
21 0 : void const * FD_RESTRICT src,
22 0 : ulong sz );
23 0 : uchar out[ 32 ];
24 0 : char buf[ 65 ] = { 0 };
25 0 : fd_f25519_tobytes( out, a );
26 0 : fd_hex_encode( buf, out, 32UL );
27 0 : FD_LOG_WARNING(( "%s: %s", name, buf ));
28 0 : FD_LOG_HEXDUMP_WARNING(( name, a, sizeof(fd_f25519_t) ));
29 0 : }
30 :
31 : /* fd_f25519_pow22523 computes r = a^(2^252-3), and returns r. */
32 : fd_f25519_t *
33 : fd_f25519_pow22523( fd_f25519_t * r,
34 2629289 : fd_f25519_t const * a ) {
35 2629289 : fd_f25519_t t0[1];
36 2629289 : fd_f25519_t t1[1];
37 2629289 : fd_f25519_t t2[1];
38 :
39 2629289 : fd_f25519_sqr( t0, a );
40 2629289 : fd_f25519_sqr( t1, t0 );
41 5258578 : for( int i=1; i< 2; i++ ) fd_f25519_sqr( t1, t1 );
42 :
43 2629289 : fd_f25519_mul( t1, a, t1 );
44 2629289 : fd_f25519_mul( t0, t0, t1 );
45 2629289 : fd_f25519_sqr( t0, t0 );
46 2629289 : fd_f25519_mul( t0, t1, t0 );
47 2629289 : fd_f25519_sqr( t1, t0 );
48 13146445 : for( int i=1; i< 5; i++ ) fd_f25519_sqr( t1, t1 );
49 :
50 2629289 : fd_f25519_mul( t0, t1, t0 );
51 2629289 : fd_f25519_sqr( t1, t0 );
52 26292890 : for( int i=1; i< 10; i++ ) fd_f25519_sqr( t1, t1 );
53 :
54 2629289 : fd_f25519_mul( t1, t1, t0 );
55 2629289 : fd_f25519_sqr( t2, t1 );
56 52585780 : for( int i=1; i< 20; i++ ) fd_f25519_sqr( t2, t2 );
57 :
58 2629289 : fd_f25519_mul( t1, t2, t1 );
59 2629289 : fd_f25519_sqr( t1, t1 );
60 26292890 : for( int i=1; i< 10; i++ ) fd_f25519_sqr( t1, t1 );
61 :
62 2629289 : fd_f25519_mul( t0, t1, t0 );
63 2629289 : fd_f25519_sqr( t1, t0 );
64 131464450 : for( int i=1; i< 50; i++ ) fd_f25519_sqr( t1, t1 );
65 :
66 2629289 : fd_f25519_mul( t1, t1, t0 );
67 2629289 : fd_f25519_sqr( t2, t1 );
68 262928900 : for( int i=1; i<100; i++ ) fd_f25519_sqr( t2, t2 );
69 :
70 2629289 : fd_f25519_mul( t1, t2, t1 );
71 2629289 : fd_f25519_sqr( t1, t1 );
72 131464450 : for( int i=1; i< 50; i++ ) fd_f25519_sqr( t1, t1 );
73 :
74 2629289 : fd_f25519_mul( t0, t1, t0 );
75 2629289 : fd_f25519_sqr( t0, t0 );
76 5258578 : for( int i=1; i< 2; i++ ) fd_f25519_sqr( t0, t0 );
77 :
78 2629289 : fd_f25519_mul(r, t0, a );
79 2629289 : return r;
80 2629289 : }
81 :
82 : #ifdef FD_HAS_S2NBIGNUM
83 :
84 : /* fd_f25519_inv computes r = 1/a, and returns r. */
85 : fd_f25519_t *
86 : fd_f25519_inv( fd_f25519_t * r,
87 260309 : fd_f25519_t const * a ) {
88 260309 : ulong x[ 4 ];
89 260309 : ulong z[ 4 ];
90 260309 : fd_f25519_tobytes( (uchar *)x, a );
91 260309 : bignum_inv_p25519( z, x );
92 260309 : fd_msan_unpoison( z, 32UL );
93 260309 : return fd_f25519_frombytes( r, (uchar const *)z );
94 260309 : }
95 :
96 : #else
97 :
98 : /* fd_f25519_inv computes r = 1/a, and returns r. */
99 : fd_f25519_t *
100 : fd_f25519_inv( fd_f25519_t * r,
101 : fd_f25519_t const * a ) {
102 : fd_f25519_t t0[1];
103 : fd_f25519_t t1[1];
104 : fd_f25519_t t2[1];
105 : fd_f25519_t t3[1];
106 :
107 : /* Compute z**-1 = z**(2**255 - 19 - 2) with the exponent as
108 : 2**255 - 21 = (2**5) * (2**250 - 1) + 11. */
109 :
110 : fd_f25519_sqr( t0, a ); /* t0 = z**2 */
111 : fd_f25519_sqr( t1, t0 );
112 : fd_f25519_sqr( t1, t1 ); /* t1 = t0**(2**2) = z**8 */
113 : fd_f25519_mul( t1, a, t1 ); /* t1 = z * t1 = z**9 */
114 : fd_f25519_mul( t0, t0, t1 ); /* t0 = t0 * t1 = z**11 -- stash t0 away for the end. */
115 : fd_f25519_sqr( t2, t0 ); /* t2 = t0**2 = z**22 */
116 : fd_f25519_mul( t1, t1, t2 ); /* t1 = t1 * t2 = z**(2**5 - 1) */
117 : fd_f25519_sqr( t2, t1 );
118 : for( int i=1; i< 5; i++ ) fd_f25519_sqr( t2, t2 ); /* t2 = t1**(2**5) = z**((2**5) * (2**5 - 1)) */
119 : fd_f25519_mul( t1, t2, t1 ); /* t1 = t1 * t2 = z**((2**5 + 1) * (2**5 - 1)) = z**(2**10 - 1) */
120 : fd_f25519_sqr( t2, t1 );
121 : for( int i=1; i< 10; i++ ) fd_f25519_sqr( t2, t2 );
122 : fd_f25519_mul( t2, t2, t1 ); /* t2 = z**(2**20 - 1) */
123 : fd_f25519_sqr( t3, t2 );
124 : for( int i=1; i< 20; i++ ) fd_f25519_sqr( t3, t3 );
125 : fd_f25519_mul( t2, t3, t2 ); /* t2 = z**(2**40 - 1) */
126 : for( int i=0; i< 10; i++ ) fd_f25519_sqr( t2, t2 ); /* t2 = z**(2**10) * (2**40 - 1) */
127 : fd_f25519_mul( t1, t2, t1 ); /* t1 = z**(2**50 - 1) */
128 : fd_f25519_sqr( t2, t1 );
129 : for( int i=1; i< 50; i++ ) fd_f25519_sqr( t2, t2 );
130 : fd_f25519_mul( t2, t2, t1 ); /* t2 = z**(2**100 - 1) */
131 : fd_f25519_sqr( t3, t2 );
132 : for( int i=1; i<100; i++ ) fd_f25519_sqr( t3, t3 );
133 : fd_f25519_mul( t2, t3, t2 ); /* t2 = z**(2**200 - 1) */
134 : fd_f25519_sqr( t2, t2 );
135 : for( int i=1; i< 50; i++ ) fd_f25519_sqr( t2, t2 ); /* t2 = z**((2**50) * (2**200 - 1) */
136 : fd_f25519_mul( t1, t2, t1 ); /* t1 = z**(2**250 - 1) */
137 : fd_f25519_sqr( t1, t1 );
138 : for( int i=1; i< 5; i++ ) fd_f25519_sqr( t1, t1 ); /* t1 = z**((2**5) * (2**250 - 1)) */
139 : return fd_f25519_mul( r, t1, t0 ); /* Recall t0 = z**11; out = z**(2**255 - 21) */
140 : }
141 :
142 : #endif /* FD_HAS_S2NBIGNUM */
143 :
144 : /* Variable time Jacobi symbol mod p, a port of
145 : secp256k1_jacobi64_maybe_var from libsecp256k1 v0.7.1
146 : (src/modinv64_impl.h, MIT license, see NOTICE). It runs Bernstein-
147 : Yang "posdivsteps" in batches of 62 on numbers stored as 5 signed
148 : 62-bit limbs, tracking the sign of the Jacobi symbol (g|f) from f mod
149 : 8 and g mod 4. See doc/safegcd_implementation.md in libsecp256k1. */
150 :
151 34870512 : #define FD_F25519_M62 (ULONG_MAX>>2)
152 :
153 : /* fd_f25519_posdivsteps_62_var runs 62 posdivsteps on f,g given their
154 : bottom 64 bits f0,g0 and returns the new eta. t is set to the
155 : transition matrix [u,v,q,r] scaled by 2^62. Bit 0 of *jac is flipped
156 : iff the Jacobi symbol (g|f) changes sign. */
157 :
158 : static long
159 : fd_f25519_posdivsteps_62_var( long eta,
160 : ulong f0,
161 : ulong g0,
162 : long t[ 4 ],
163 9267253 : int * jac ) {
164 9267253 : ulong u = 1UL, v = 0UL, q = 0UL, r = 1UL;
165 9267253 : ulong f = f0, g = g0;
166 9267253 : int i = 62;
167 9267253 : int j = *jac;
168 :
169 164279256 : for(;;) {
170 : /* Divide g by 2 as often as possible (at most i times). When
171 : dividing by an odd power of 2, the symbol flips if f mod 8 is 3
172 : or 5. */
173 164279256 : int zeros = fd_ulong_find_lsb( g | (ULONG_MAX << i) );
174 164279256 : g >>= zeros;
175 164279256 : u <<= zeros;
176 164279256 : v <<= zeros;
177 164279256 : eta -= zeros;
178 164279256 : i -= zeros;
179 164279256 : j ^= (int)( (ulong)zeros & ((f>>1) ^ (f>>2)) );
180 164279256 : if( !i ) break;
181 :
182 155012003 : ulong m;
183 155012003 : ulong w;
184 155012003 : int limit;
185 155012003 : if( eta<0L ) {
186 : /* Swap f,g. The symbol flips if both are 3 mod 4. Then cancel
187 : up to 6 bottom bits of g with a multiple of f. */
188 149837612 : ulong tmp;
189 149837612 : eta = -eta;
190 149837612 : tmp = f; f = g; g = tmp;
191 149837612 : tmp = u; u = q; q = tmp;
192 149837612 : tmp = v; v = r; r = tmp;
193 149837612 : j ^= (int)( (f & g)>>1 );
194 149837612 : limit = fd_int_min( (int)eta+1, i );
195 149837612 : m = (ULONG_MAX >> (64-limit)) & 63UL;
196 149837612 : w = (f*g*(f*f-2UL)) & m;
197 149837612 : } else {
198 : /* Cancel up to 4 bottom bits of g with a multiple of f. */
199 5174391 : limit = fd_int_min( (int)eta+1, i );
200 5174391 : m = (ULONG_MAX >> (64-limit)) & 15UL;
201 5174391 : w = f + (((f+1UL) & 4UL)<<1);
202 5174391 : w = (-w*g) & m;
203 5174391 : }
204 155012003 : g += f*w;
205 155012003 : q += u*w;
206 155012003 : r += v*w;
207 155012003 : }
208 :
209 9267253 : t[0] = (long)u; t[1] = (long)v; t[2] = (long)q; t[3] = (long)r;
210 9267253 : *jac = j;
211 9267253 : return eta;
212 9267253 : }
213 :
214 : /* fd_f25519_update_fg_62_var computes [f,g] = t [f,g] / 2^62 on the
215 : bottom len limbs of f and g. */
216 :
217 : static void
218 : fd_f25519_update_fg_62_var( int len,
219 : long * f,
220 : long * g,
221 9267253 : long const t[ 4 ] ) {
222 9267253 : long const u = t[0], v = t[1], q = t[2], r = t[3];
223 9267253 : int128 cf = (int128)u*(int128)f[0] + (int128)v*(int128)g[0];
224 9267253 : int128 cg = (int128)q*(int128)f[0] + (int128)r*(int128)g[0];
225 9267253 : cf >>= 62;
226 9267253 : cg >>= 62;
227 25244929 : for( int i=1; i<len; i++ ) {
228 15977676 : cf += (int128)u*(int128)f[i] + (int128)v*(int128)g[i];
229 15977676 : cg += (int128)q*(int128)f[i] + (int128)r*(int128)g[i];
230 15977676 : f[i-1] = (long)( (ulong)cf & FD_F25519_M62 ); cf >>= 62;
231 15977676 : g[i-1] = (long)( (ulong)cg & FD_F25519_M62 ); cg >>= 62;
232 15977676 : }
233 9267253 : f[len-1] = (long)cf;
234 9267253 : g[len-1] = (long)cg;
235 9267253 : }
236 :
237 : /* fd_f25519_jacobi_var returns the Jacobi symbol (x|p) in {-1,1} for x
238 : in [1,p) given as 5 non-negative 62-bit limbs, or 0 if it did not
239 : converge within 25*62 posdivsteps (never observed for random input,
240 : typically ~12 batches are needed). */
241 :
242 : static int
243 728790 : fd_f25519_jacobi_var( long const x[ 5 ] ) {
244 728790 : long f[5] = { (1L<<62)-19L, (1L<<62)-1L, (1L<<62)-1L, (1L<<62)-1L, 127L }; /* p */
245 728790 : long g[5] = { x[0], x[1], x[2], x[3], x[4] };
246 728790 : int len = 5;
247 728790 : long eta = -1L;
248 728790 : int jac = 0;
249 :
250 9267382 : for( int cnt=0; cnt<25; cnt++ ) {
251 9267253 : long t[4];
252 9267253 : eta = fd_f25519_posdivsteps_62_var( eta, (ulong)f[0] | ((ulong)f[1]<<62), (ulong)g[0] | ((ulong)g[1]<<62), t, &jac );
253 9267253 : fd_f25519_update_fg_62_var( len, f, g, t );
254 :
255 : /* Done when f==1 */
256 9267253 : if( f[0]==1L ) {
257 749592 : long hi = 0L;
258 821586 : for( int j=1; j<len; j++ ) hi |= f[j];
259 749592 : if( !hi ) return 1 - 2*(jac & 1);
260 749592 : }
261 :
262 : /* Drop the top limb when it is zero in both f and g */
263 8538592 : if( len>1 && !(f[len-1] | g[len-1]) ) len--;
264 8538592 : }
265 129 : return 0;
266 728790 : }
267 :
268 : int
269 728832 : fd_f25519_is_square_var( fd_f25519_t const * a ) {
270 728832 : uchar buf[32];
271 728832 : fd_f25519_tobytes( buf, a ); /* canonical, in [0,p) */
272 728832 : ulong a0 = fd_ulong_load_8_fast( buf );
273 728832 : ulong a1 = fd_ulong_load_8_fast( buf+ 8 );
274 728832 : ulong a2 = fd_ulong_load_8_fast( buf+16 );
275 728832 : ulong a3 = fd_ulong_load_8_fast( buf+24 );
276 728832 : if( FD_UNLIKELY( !(a0|a1|a2|a3) ) ) return 1;
277 :
278 728790 : long x[5] = { (long)( a0 & FD_F25519_M62 ),
279 728790 : (long)( ((a0>>62) | (a1<< 2)) & FD_F25519_M62 ),
280 728790 : (long)( ((a1>>60) | (a2<< 4)) & FD_F25519_M62 ),
281 728790 : (long)( ((a2>>58) | (a3<< 6)) & FD_F25519_M62 ),
282 728790 : (long)( a3>>56 ) };
283 728790 : int jac = fd_f25519_jacobi_var( x );
284 728790 : if( FD_UNLIKELY( !jac ) ) {
285 : /* Did not converge: a is a nonzero square iff sqrt(a/1) exists */
286 129 : fd_f25519_t r[1];
287 129 : return fd_f25519_sqrt_ratio( r, a, fd_f25519_one );
288 129 : }
289 728661 : return jac>0;
290 728790 : }
291 :
292 : #undef FD_F25519_M62
293 :
294 : /* fd_f25519_sqrt_ratio computes r = (u * v^3) * (u * v^7)^((p-5)/8),
295 : returns 0 on success, 1 on failure. */
296 : int
297 : fd_f25519_sqrt_ratio( fd_f25519_t * r,
298 : fd_f25519_t const * u,
299 2515550 : fd_f25519_t const * v ) {
300 : /* r = (u * v^3) * (u * v^7)^((p-5)/8) */
301 2515550 : fd_f25519_t v2[1]; fd_f25519_sqr( v2, v );
302 2515550 : fd_f25519_t v3[1]; fd_f25519_mul( v3, v2, v );
303 2515550 : fd_f25519_t uv3[1]; fd_f25519_mul( uv3, u, v3 );
304 2515550 : fd_f25519_t v6[1]; fd_f25519_sqr( v6, v3 );
305 2515550 : fd_f25519_t v7[1]; fd_f25519_mul( v7, v6, v );
306 2515550 : fd_f25519_t uv7[1]; fd_f25519_mul( uv7, u, v7 );
307 2515550 : fd_f25519_pow22523( r, uv7 );
308 2515550 : fd_f25519_mul ( r, r, uv3 );
309 :
310 : /* check = v * r^2 */
311 2515550 : fd_f25519_t check[1];
312 2515550 : fd_f25519_sqr( check, r );
313 2515550 : fd_f25519_mul( check, check, v );
314 :
315 : /* (correct_sign_sqrt) check == u
316 : (flipped_sign_sqrt) check == !u
317 : (flipped_sign_sqrt_i) check == (!u * SQRT_M1) */
318 2515550 : fd_f25519_t u_neg[1]; fd_f25519_neg( u_neg, u );
319 2515550 : fd_f25519_t u_neg_sqrtm1[1]; fd_f25519_mul( u_neg_sqrtm1, u_neg, fd_f25519_sqrtm1 );
320 2515550 : int correct_sign_sqrt = fd_f25519_eq( check, u );
321 2515550 : int flipped_sign_sqrt = fd_f25519_eq( check, u_neg );
322 2515550 : int flipped_sign_sqrt_i = fd_f25519_eq( check, u_neg_sqrtm1 );
323 :
324 : /* r_prime = SQRT_M1 * r */
325 2515550 : fd_f25519_t r_prime[1];
326 2515550 : fd_f25519_mul( r_prime, r, fd_f25519_sqrtm1 );
327 :
328 : /* r = CT_SELECT(r_prime IF flipped_sign_sqrt | flipped_sign_sqrt_i ELSE r) */
329 2515550 : fd_f25519_if( r, flipped_sign_sqrt|flipped_sign_sqrt_i, r_prime, r );
330 2515550 : fd_f25519_abs( r, r );
331 2515550 : return correct_sign_sqrt|flipped_sign_sqrt;
332 2515550 : }
|
