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