Line data Source code
1 : #include "fd_big_mod_exp.h"
2 :
3 : #if FD_HAS_S2NBIGNUM
4 :
5 : #include <stdint.h>
6 : #include <s2n-bignum.h>
7 :
8 : /* s2n-bignum has ADX variants for various sizes, and generic non-ADX
9 : bignum_amontmul, bignum_amontsqr.
10 : ADX code is ~2x faster than non-ADX, so worth the extra complexity. */
11 : #if defined(__ADX__)
12 9600 : #define FD_BIG_MOD_EXP_ADX 1
13 : #else
14 19200 : #define FD_BIG_MOD_EXP_ADX 0
15 : #endif
16 :
17 : /* Max operand size in 64-bit words (512 bytes / 8). */
18 : #define FD_BIG_MOD_EXP_MAX_WORDS (FD_BIG_MOD_EXP_MAX_BYTES/8UL)
19 :
20 : /* Largest exponent window, in bits. */
21 :
22 57600 : #define FD_BIG_MOD_EXP_MAX_WINDOW (4UL)
23 : #define FD_BIG_MOD_EXP_TBL_CNT (1UL<<FD_BIG_MOD_EXP_MAX_WINDOW)
24 :
25 : /* largest k{sqr,mul} temp (kmul_32_64) */
26 : #define FD_BIG_MOD_EXP_KT_WORDS (96UL)
27 :
28 : int
29 : fd_big_mod_exp_modulus_is_valid( uchar const * mod,
30 69339 : ulong mod_len ) {
31 69339 : if( FD_UNLIKELY( mod_len==0UL ) ) return 0; /* zero */
32 69330 : if( FD_UNLIKELY( !(mod[0] & 1U) ) ) return 0; /* even (includes zero) */
33 : /* odd: the only remaining value <= 1 is exactly 1 */
34 69279 : if( FD_UNLIKELY( mod[0]==1U ) ) {
35 429 : int is_one = 1;
36 471 : for( ulong i=1UL; i<mod_len; i++ ) if( mod[i] ) { is_one = 0; break; }
37 429 : if( is_one ) return 0; /* one */
38 429 : }
39 69249 : return 1;
40 69279 : }
41 :
42 : FD_FN_PURE static inline ulong
43 : fd_big_mod_exp_bits( ulong const * x,
44 : ulong pos,
45 400968 : ulong n ) {
46 400968 : ulong i = pos>>6;
47 400968 : ulong sh = pos & 63UL;
48 400968 : ulong v = x[i]>>sh;
49 400968 : if( sh+n>64UL ) v |= x[i+1]<<(64UL-sh);
50 400968 : return v & ((1UL<<n)-1UL);
51 400968 : }
52 :
53 : /* Montgomery domain parameters plus the fast path's scratch. */
54 :
55 : struct fd_big_mod_exp_mont {
56 : ulong const * m; /* modulus, k words */
57 : ulong k;
58 : ulong w; /* -1/m[0] mod 2^64 */
59 : int fast;
60 : ulong * prod; /* 2k words */
61 : ulong * kt; /* k{sqr,mul} temporary */
62 : };
63 :
64 : typedef struct fd_big_mod_exp_mont fd_big_mod_exp_mont_t;
65 :
66 : static inline void
67 : fd_big_mod_exp_mont_sqr( fd_big_mod_exp_mont_t const * mt,
68 : ulong * out,
69 1408998 : ulong const * x ) {
70 1408998 : ulong k = mt->k;
71 469666 : #if FD_BIG_MOD_EXP_ADX
72 469666 : if( FD_LIKELY( mt->fast ) ) {
73 20778 : if ( k== 8UL ) bignum_sqr_8_16 ( mt->prod, x );
74 18457 : else if( k==16UL ) bignum_ksqr_16_32( mt->prod, x, mt->kt );
75 16750 : else if( k==32UL ) bignum_ksqr_32_64( mt->prod, x, mt->kt );
76 16097 : else bignum_sqr( 2UL*k, mt->prod, k, x );
77 20778 : ulong cy = bignum_emontredc_8n( k, mt->prod, mt->m, mt->w );
78 20778 : bignum_optsub( k, out, mt->prod+k, cy, mt->m );
79 20778 : return;
80 20778 : }
81 448888 : #endif
82 1388220 : bignum_amontsqr( k, out, x, mt->m );
83 1388220 : }
84 :
85 : static inline void
86 : fd_big_mod_exp_mont_mul( fd_big_mod_exp_mont_t const * mt,
87 : ulong * out,
88 : ulong const * x,
89 410274 : ulong const * y ) {
90 410274 : ulong k = mt->k;
91 136758 : #if FD_BIG_MOD_EXP_ADX
92 136758 : if( FD_LIKELY( mt->fast ) ) {
93 4486 : if ( k== 8UL ) bignum_mul_8_16 ( mt->prod, x, y );
94 3865 : else if( k==16UL ) bignum_kmul_16_32( mt->prod, x, y, mt->kt );
95 3449 : else if( k==32UL ) bignum_kmul_32_64( mt->prod, x, y, mt->kt );
96 3265 : else bignum_mul( 2UL*k, mt->prod, k, x, k, y );
97 4486 : ulong cy = bignum_emontredc_8n( k, mt->prod, mt->m, mt->w );
98 4486 : bignum_optsub( k, out, mt->prod+k, cy, mt->m );
99 4486 : return;
100 4486 : }
101 132272 : #endif
102 405788 : bignum_amontmul( k, out, x, y, mt->m );
103 405788 : }
104 :
105 : /* r = a mod m, all k words, a any value, m>1.
106 : Specialized implementation to keep a low CU cost for the syscall. */
107 : static void
108 : fd_big_mod_exp_reduce( ulong * r,
109 : ulong const * a,
110 : ulong const * m,
111 21696 : ulong k ) {
112 403890 : ulong km = k; while( !m[km-1UL] ) km--; /* m>1, so km>=1 */
113 463404 : ulong ka = k; while( ka && !a[ka-1UL] ) ka--;
114 :
115 21696 : fd_memset( r, 0, k*sizeof(ulong) );
116 :
117 : /* a < m: nothing to reduce */
118 21696 : if( ka<km ) { fd_memcpy( r, a, ka*sizeof(ulong) ); return; }
119 12093 : if( ka==km ) {
120 3483 : ulong i = km;
121 22947 : while( i && a[i-1UL]==m[i-1UL] ) i--;
122 3483 : if( i && a[i-1UL]<m[i-1UL] ) { fd_memcpy( r, a, ka*sizeof(ulong) ); return; }
123 3483 : }
124 :
125 10824 : ulong mn[ FD_BIG_MOD_EXP_MAX_WORDS ]; /* m << s, km words */
126 10824 : ulong an[ FD_BIG_MOD_EXP_MAX_WORDS+1UL ]; /* a << s, ka+1 words */
127 :
128 10824 : ulong s = 63UL - (ulong)fd_ulong_find_msb( m[km-1UL] );
129 10824 : if( s ) {
130 53016 : for( ulong i=km; i-->1UL; ) mn[i] = (m[i]<<s) | (m[i-1UL]>>(64UL-s));
131 5643 : mn[0] = m[0]<<s;
132 5643 : an[ka] = a[ka-1UL]>>(64UL-s);
133 114540 : for( ulong i=ka; i-->1UL; ) an[i] = (a[i]<<s) | (a[i-1UL]>>(64UL-s));
134 5643 : an[0] = a[0]<<s;
135 5643 : } else {
136 5181 : fd_memcpy( mn, m, km*sizeof(ulong) );
137 5181 : fd_memcpy( an, a, ka*sizeof(ulong) );
138 5181 : an[ka] = 0UL;
139 5181 : }
140 :
141 10824 : if( km==1UL ) {
142 2922 : ulong d = mn[0];
143 2922 : uint128 rem = 0;
144 38592 : for( ulong j=ka+1UL; j-->0UL; ) rem = ( (rem<<64) | an[j] ) % d; /* rem<d, so no overflow */
145 2922 : r[0] = (ulong)rem>>s;
146 2922 : return;
147 2922 : }
148 :
149 7902 : ulong dh = mn[km-1UL];
150 7902 : ulong dl = mn[km-2UL];
151 130740 : for( ulong j=ka-km+1UL; j-->0UL; ) {
152 : /* quotient digit estimate for the window an[j..j+km]; an[j+km]<=dh */
153 122838 : ulong qhat;
154 122838 : uint128 rhat;
155 122838 : if( an[j+km]>=dh ) {
156 0 : qhat = ~0UL;
157 0 : rhat = (uint128)an[j+km-1UL] + dh;
158 122838 : } else {
159 122838 : uint128 num = ( (uint128)an[j+km]<<64 ) | an[j+km-1UL];
160 122838 : qhat = (ulong)( num/dh );
161 122838 : rhat = num%dh;
162 122838 : }
163 163692 : while( !(rhat>>64) && (uint128)qhat*dl > ( (rhat<<64) | an[j+km-2UL] ) ) {
164 40854 : qhat--;
165 40854 : rhat += dh;
166 40854 : }
167 :
168 : /* an[j..j+km] -= qhat*mn */
169 122838 : ulong carry = 0UL;
170 122838 : ulong borrow = 0UL;
171 1436667 : for( ulong i=0UL; i<km; i++ ) {
172 1313829 : uint128 p = (uint128)qhat*mn[i] + carry;
173 1313829 : carry = (ulong)( p>>64 );
174 1313829 : uint128 d = (uint128)an[j+i] - (ulong)p - borrow;
175 1313829 : an[j+i] = (ulong)d;
176 1313829 : borrow = (ulong)( d>>64 ) & 1UL;
177 1313829 : }
178 122838 : uint128 d = (uint128)an[j+km] - carry - borrow;
179 122838 : an[j+km] = (ulong)d;
180 :
181 : /* estimate one too large: add mn back */
182 122838 : if( d>>64 ) {
183 594 : ulong c = 0UL;
184 16821 : for( ulong i=0UL; i<km; i++ ) {
185 16227 : uint128 t = (uint128)an[j+i] + mn[i] + c;
186 16227 : an[j+i] = (ulong)t;
187 16227 : c = (ulong)( t>>64 );
188 16227 : }
189 594 : an[j+km] += c;
190 594 : }
191 122838 : }
192 :
193 : /* remainder is an[0..km), with an[km]==0 */
194 7902 : if( s ) {
195 55047 : for( ulong i=0UL; i<km; i++ ) r[i] = (an[i]>>s) | (an[i+1UL]<<(64UL-s));
196 4065 : } else {
197 4065 : fd_memcpy( r, an, km*sizeof(ulong) );
198 4065 : }
199 7902 : }
200 :
201 : int
202 : fd_big_mod_exp( uchar * out,
203 : uchar const * base, ulong base_len,
204 : uchar const * exp, ulong exp_len,
205 43065 : uchar const * mod, ulong mod_len ) {
206 :
207 43065 : if( FD_UNLIKELY( !fd_big_mod_exp_modulus_is_valid( mod, mod_len ) ) )
208 33 : return FD_BIG_MOD_EXP_ERR_MODULUS;
209 :
210 43032 : ulong k = ( fd_ulong_max( base_len, mod_len ) + 7UL )/8UL; /* [1,MAX_WORDS] */
211 43032 : ulong e = ( exp_len + 7UL )/8UL; /* [0,MAX_WORDS] */
212 :
213 43032 : ulong a [ FD_BIG_MOD_EXP_MAX_WORDS ]; /* base, zero extended to k */
214 43032 : ulong p [ FD_BIG_MOD_EXP_MAX_WORDS ]; /* exponent, zero extended to e */
215 43032 : ulong m [ FD_BIG_MOD_EXP_MAX_WORDS ]; /* modulus, zero extended to k */
216 43032 : ulong r2 [ FD_BIG_MOD_EXP_MAX_WORDS ]; /* 2^{128k} mod m */
217 43032 : ulong u [ FD_BIG_MOD_EXP_MAX_WORDS ]; /* accumulator, ping */
218 43032 : ulong v [ FD_BIG_MOD_EXP_MAX_WORDS ]; /* accumulator, pong */
219 43032 : ulong tbl[ FD_BIG_MOD_EXP_TBL_CNT ][ FD_BIG_MOD_EXP_MAX_WORDS ]; /* mont(base^i) */
220 43032 : ulong prd[ 2UL*FD_BIG_MOD_EXP_MAX_WORDS ]; /* fast path: product */
221 43032 : ulong kt [ FD_BIG_MOD_EXP_KT_WORDS ]; /* fast path: k{sqr,mul} temp */
222 :
223 43032 : fd_memset( a, 0, k*sizeof(ulong) );
224 43032 : fd_memset( m, 0, k*sizeof(ulong) );
225 43032 : fd_memcpy( (uchar *)a, base, base_len );
226 43032 : fd_memcpy( (uchar *)m, mod, mod_len );
227 43032 : if( FD_LIKELY( e ) ) {
228 40833 : fd_memset( p, 0, e*sizeof(ulong) );
229 40833 : fd_memcpy( (uchar *)p, exp, exp_len );
230 40833 : }
231 :
232 43032 : ulong nbits = e ? bignum_bitsize( e, p ) : 0UL;
233 :
234 : /* base^0 == 1 for every base, 0 included. m>1, so 1 mod m == 1. */
235 43032 : if( FD_UNLIKELY( !nbits ) ) {
236 6936 : fd_memset( out, 0, mod_len );
237 6936 : out[0] = (uchar)1;
238 6936 : return FD_BIG_MOD_EXP_SUCCESS;
239 6936 : }
240 :
241 : /* base^1 == base mod m: one reduction, no Montgomery setup. */
242 36096 : if( FD_UNLIKELY( nbits==1UL ) ) {
243 21696 : fd_big_mod_exp_reduce( u, a, m, k );
244 21696 : fd_memcpy( out, (uchar const *)u, mod_len );
245 21696 : return FD_BIG_MOD_EXP_SUCCESS;
246 21696 : }
247 :
248 : /* Fixed window. */
249 14400 : ulong w = 1UL;
250 14400 : ulong best = nbits + nbits; /* cost at w==1 */
251 57600 : for( ulong t=2UL; t<=FD_BIG_MOD_EXP_MAX_WINDOW; t++ ) {
252 43200 : ulong cost = ((1UL<<t)-2UL) + nbits + (nbits+t-1UL)/t;
253 43200 : if( cost<best ) { best = cost; w = t; }
254 43200 : }
255 :
256 14400 : fd_big_mod_exp_mont_t mt[1];
257 14400 : mt->m = m;
258 14400 : mt->k = k;
259 14400 : mt->prod = prd;
260 14400 : mt->kt = kt;
261 14400 : mt->fast = FD_BIG_MOD_EXP_ADX && k>=8UL && !(k & 7UL);
262 14400 : mt->w = mt->fast ? word_negmodinv( m[0] ) : 0UL;
263 :
264 14400 : bignum_amontifier( k, r2, m, v /* scratch, >=k words */ );
265 14400 : fd_big_mod_exp_mont_mul( mt, tbl[1], r2, a );
266 125766 : for( ulong i=2UL; i<(1UL<<w); i++ ) {
267 111366 : fd_big_mod_exp_mont_mul( mt, tbl[i], tbl[i-1], tbl[1] );
268 111366 : }
269 :
270 14400 : ulong top = nbits % w;
271 14400 : if( !top ) top = w;
272 14400 : ulong pos = nbits - top; /* a multiple of w */
273 :
274 14400 : ulong * cur = u;
275 14400 : ulong * nxt = v;
276 14400 : fd_memcpy( cur, tbl[ fd_big_mod_exp_bits( p, pos, top ) ], k*sizeof(ulong) );
277 :
278 400968 : while( pos ) {
279 386568 : pos -= w;
280 1795566 : for( ulong j=0UL; j<w; j++ ) {
281 1408998 : fd_big_mod_exp_mont_sqr( mt, nxt, cur );
282 1408998 : ulong * s = cur; cur = nxt; nxt = s;
283 1408998 : }
284 386568 : ulong d = fd_big_mod_exp_bits( p, pos, w );
285 386568 : if( d ) {
286 284508 : fd_big_mod_exp_mont_mul( mt, nxt, cur, tbl[d] );
287 284508 : ulong * s = cur; cur = nxt; nxt = s;
288 284508 : }
289 386568 : }
290 :
291 14400 : bignum_demont( k, nxt, cur, m );
292 14400 : fd_memcpy( out, (uchar const *)nxt, mod_len );
293 :
294 14400 : return FD_BIG_MOD_EXP_SUCCESS;
295 36096 : }
296 :
297 : #endif /* FD_HAS_S2NBIGNUM */
|
