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 28951 : fd_f25519_t const * a ) {
35 28951 : fd_f25519_t t0[1];
36 28951 : fd_f25519_t t1[1];
37 28951 : fd_f25519_t t2[1];
38 :
39 28951 : fd_f25519_sqr( t0, a );
40 28951 : fd_f25519_sqr( t1, t0 );
41 57888 : for( int i=1; i< 2; i++ ) fd_f25519_sqr( t1, t1 );
42 :
43 28951 : fd_f25519_mul( t1, a, t1 );
44 28951 : fd_f25519_mul( t0, t0, t1 );
45 28951 : fd_f25519_sqr( t0, t0 );
46 28951 : fd_f25519_mul( t0, t1, t0 );
47 28951 : fd_f25519_sqr( t1, t0 );
48 144827 : for( int i=1; i< 5; i++ ) fd_f25519_sqr( t1, t1 );
49 :
50 28951 : fd_f25519_mul( t0, t1, t0 );
51 28951 : fd_f25519_sqr( t1, t0 );
52 286459 : for( int i=1; i< 10; i++ ) fd_f25519_sqr( t1, t1 );
53 :
54 28951 : fd_f25519_mul( t1, t1, t0 );
55 28951 : fd_f25519_sqr( t2, t1 );
56 562319 : for( int i=1; i< 20; i++ ) fd_f25519_sqr( t2, t2 );
57 :
58 28951 : fd_f25519_mul( t1, t2, t1 );
59 28951 : fd_f25519_sqr( t1, t1 );
60 286567 : for( int i=1; i< 10; i++ ) fd_f25519_sqr( t1, t1 );
61 :
62 28951 : fd_f25519_mul( t0, t1, t0 );
63 28951 : fd_f25519_sqr( t1, t0 );
64 1342445 : for( int i=1; i< 50; i++ ) fd_f25519_sqr( t1, t1 );
65 :
66 28951 : fd_f25519_mul( t1, t1, t0 );
67 28951 : fd_f25519_sqr( t2, t1 );
68 2518801 : for( int i=1; i<100; i++ ) fd_f25519_sqr( t2, t2 );
69 :
70 28951 : fd_f25519_mul( t1, t2, t1 );
71 28951 : fd_f25519_sqr( t1, t1 );
72 1340583 : for( int i=1; i< 50; i++ ) fd_f25519_sqr( t1, t1 );
73 :
74 28951 : fd_f25519_mul( t0, t1, t0 );
75 28951 : fd_f25519_sqr( t0, t0 );
76 57974 : for( int i=1; i< 2; i++ ) fd_f25519_sqr( t0, t0 );
77 :
78 28951 : fd_f25519_mul(r, t0, a );
79 28951 : return r;
80 28951 : }
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 283 : fd_f25519_t const * a ) {
88 283 : ulong x[ 4 ];
89 283 : ulong z[ 4 ];
90 283 : fd_f25519_tobytes( (uchar *)x, a );
91 283 : bignum_inv_p25519( z, x );
92 283 : fd_msan_unpoison( z, 32UL );
93 283 : return fd_f25519_frombytes( r, (uchar const *)z );
94 283 : }
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 93928 : #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 25048 : int * jac ) {
164 25048 : ulong u = 1UL, v = 0UL, q = 0UL, r = 1UL;
165 25048 : ulong f = f0, g = g0;
166 25048 : int i = 62;
167 25048 : int j = *jac;
168 :
169 440503 : 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 440503 : int zeros = fd_ulong_find_lsb( g | (ULONG_MAX << i) );
174 440503 : g >>= zeros;
175 440503 : u <<= zeros;
176 440503 : v <<= zeros;
177 440503 : eta -= zeros;
178 440503 : i -= zeros;
179 440503 : j ^= (int)( (ulong)zeros & ((f>>1) ^ (f>>2)) );
180 440503 : if( !i ) break;
181 :
182 415455 : ulong m;
183 415455 : ulong w;
184 415455 : int limit;
185 415455 : 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 402329 : ulong tmp;
189 402329 : eta = -eta;
190 402329 : tmp = f; f = g; g = tmp;
191 402329 : tmp = u; u = q; q = tmp;
192 402329 : tmp = v; v = r; r = tmp;
193 402329 : j ^= (int)( (f & g)>>1 );
194 402329 : limit = fd_int_min( (int)eta+1, i );
195 402329 : m = (ULONG_MAX >> (64-limit)) & 63UL;
196 402329 : w = (f*g*(f*f-2UL)) & m;
197 402329 : } else {
198 : /* Cancel up to 4 bottom bits of g with a multiple of f. */
199 13126 : limit = fd_int_min( (int)eta+1, i );
200 13126 : m = (ULONG_MAX >> (64-limit)) & 15UL;
201 13126 : w = f + (((f+1UL) & 4UL)<<1);
202 13126 : w = (-w*g) & m;
203 13126 : }
204 415455 : g += f*w;
205 415455 : q += u*w;
206 415455 : r += v*w;
207 415455 : }
208 :
209 25048 : t[0] = (long)u; t[1] = (long)v; t[2] = (long)q; t[3] = (long)r;
210 25048 : *jac = j;
211 25048 : return eta;
212 25048 : }
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 25048 : long const t[ 4 ] ) {
222 25048 : long const u = t[0], v = t[1], q = t[2], r = t[3];
223 25048 : int128 cf = (int128)u*(int128)f[0] + (int128)v*(int128)g[0];
224 25048 : int128 cg = (int128)q*(int128)f[0] + (int128)r*(int128)g[0];
225 25048 : cf >>= 62;
226 25048 : cg >>= 62;
227 68010 : for( int i=1; i<len; i++ ) {
228 42962 : cf += (int128)u*(int128)f[i] + (int128)v*(int128)g[i];
229 42962 : cg += (int128)q*(int128)f[i] + (int128)r*(int128)g[i];
230 42962 : f[i-1] = (long)( (ulong)cf & FD_F25519_M62 ); cf >>= 62;
231 42962 : g[i-1] = (long)( (ulong)cg & FD_F25519_M62 ); cg >>= 62;
232 42962 : }
233 25048 : f[len-1] = (long)cf;
234 25048 : g[len-1] = (long)cg;
235 25048 : }
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 2001 : fd_f25519_jacobi_var( long const x[ 5 ] ) {
244 2001 : long f[5] = { (1L<<62)-19L, (1L<<62)-1L, (1L<<62)-1L, (1L<<62)-1L, 127L }; /* p */
245 2001 : long g[5] = { x[0], x[1], x[2], x[3], x[4] };
246 2001 : int len = 5;
247 2001 : long eta = -1L;
248 2001 : int jac = 0;
249 :
250 25049 : for( int cnt=0; cnt<25; cnt++ ) {
251 25048 : long t[4];
252 25048 : eta = fd_f25519_posdivsteps_62_var( eta, (ulong)f[0] | ((ulong)f[1]<<62), (ulong)g[0] | ((ulong)g[1]<<62), t, &jac );
253 25048 : fd_f25519_update_fg_62_var( len, f, g, t );
254 :
255 : /* Done when f==1 */
256 25048 : if( f[0]==1L ) {
257 1987 : long hi = 0L;
258 1987 : for( int j=1; j<len; j++ ) hi |= f[j];
259 2000 : if( !hi ) return 1 - 2*(jac & 1);
260 1987 : }
261 :
262 : /* Drop the top limb when it is zero in both f and g */
263 23048 : if( len>1 && !(f[len-1] | g[len-1]) ) len--;
264 23048 : }
265 1 : return 0;
266 2001 : }
267 :
268 : int
269 2001 : fd_f25519_is_square_var( fd_f25519_t const * a ) {
270 2001 : uchar buf[32];
271 2001 : fd_f25519_tobytes( buf, a ); /* canonical, in [0,p) */
272 2001 : ulong a0 = fd_ulong_load_8_fast( buf );
273 2001 : ulong a1 = fd_ulong_load_8_fast( buf+ 8 );
274 2001 : ulong a2 = fd_ulong_load_8_fast( buf+16 );
275 2001 : ulong a3 = fd_ulong_load_8_fast( buf+24 );
276 2001 : if( FD_UNLIKELY( !(a0|a1|a2|a3) ) ) return 1;
277 :
278 2001 : long x[5] = { (long)( a0 & FD_F25519_M62 ),
279 2001 : (long)( ((a0>>62) | (a1<< 2)) & FD_F25519_M62 ),
280 2001 : (long)( ((a1>>60) | (a2<< 4)) & FD_F25519_M62 ),
281 2001 : (long)( ((a2>>58) | (a3<< 6)) & FD_F25519_M62 ),
282 2001 : (long)( a3>>56 ) };
283 2001 : int jac = fd_f25519_jacobi_var( x );
284 2001 : if( FD_UNLIKELY( !jac ) ) {
285 : /* Did not converge: a is a nonzero square iff sqrt(a/1) exists */
286 0 : fd_f25519_t r[1];
287 0 : return fd_f25519_sqrt_ratio( r, a, fd_f25519_one );
288 0 : }
289 2001 : return jac>0;
290 2001 : }
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 28882 : fd_f25519_t const * v ) {
300 : /* r = (u * v^3) * (u * v^7)^((p-5)/8) */
301 28882 : fd_f25519_t v2[1]; fd_f25519_sqr( v2, v );
302 28882 : fd_f25519_t v3[1]; fd_f25519_mul( v3, v2, v );
303 28882 : fd_f25519_t uv3[1]; fd_f25519_mul( uv3, u, v3 );
304 28882 : fd_f25519_t v6[1]; fd_f25519_sqr( v6, v3 );
305 28882 : fd_f25519_t v7[1]; fd_f25519_mul( v7, v6, v );
306 28882 : fd_f25519_t uv7[1]; fd_f25519_mul( uv7, u, v7 );
307 28882 : fd_f25519_pow22523( r, uv7 );
308 28882 : fd_f25519_mul ( r, r, uv3 );
309 :
310 : /* check = v * r^2 */
311 28882 : fd_f25519_t check[1];
312 28882 : fd_f25519_sqr( check, r );
313 28882 : 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 28882 : fd_f25519_t u_neg[1]; fd_f25519_neg( u_neg, u );
319 28882 : fd_f25519_t u_neg_sqrtm1[1]; fd_f25519_mul( u_neg_sqrtm1, u_neg, fd_f25519_sqrtm1 );
320 28882 : int correct_sign_sqrt = fd_f25519_eq( check, u );
321 28882 : int flipped_sign_sqrt = fd_f25519_eq( check, u_neg );
322 28882 : int flipped_sign_sqrt_i = fd_f25519_eq( check, u_neg_sqrtm1 );
323 :
324 : /* r_prime = SQRT_M1 * r */
325 28882 : fd_f25519_t r_prime[1];
326 28882 : 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 28882 : fd_f25519_if( r, flipped_sign_sqrt|flipped_sign_sqrt_i, r_prime, r );
330 28882 : fd_f25519_abs( r, r );
331 28882 : return correct_sign_sqrt|flipped_sign_sqrt;
332 28882 : }
|
