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 0 : #define FD_BIG_MOD_EXP_ADX 1
13 : #else
14 : #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 0 : #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 0 : ulong mod_len ) {
31 0 : if( FD_UNLIKELY( mod_len==0UL ) ) return 0; /* zero */
32 0 : if( FD_UNLIKELY( !(mod[0] & 1U) ) ) return 0; /* even (includes zero) */
33 : /* odd: the only remaining value <= 1 is exactly 1 */
34 0 : if( FD_UNLIKELY( mod[0]==1U ) ) {
35 0 : int is_one = 1;
36 0 : for( ulong i=1UL; i<mod_len; i++ ) if( mod[i] ) { is_one = 0; break; }
37 0 : if( is_one ) return 0; /* one */
38 0 : }
39 0 : return 1;
40 0 : }
41 :
42 : FD_FN_PURE static inline ulong
43 : fd_big_mod_exp_bits( ulong const * x,
44 : ulong pos,
45 0 : ulong n ) {
46 0 : ulong i = pos>>6;
47 0 : ulong sh = pos & 63UL;
48 0 : ulong v = x[i]>>sh;
49 0 : if( sh+n>64UL ) v |= x[i+1]<<(64UL-sh);
50 0 : return v & ((1UL<<n)-1UL);
51 0 : }
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 0 : ulong const * x ) {
70 0 : ulong k = mt->k;
71 0 : #if FD_BIG_MOD_EXP_ADX
72 0 : if( FD_LIKELY( mt->fast ) ) {
73 0 : if ( k== 8UL ) bignum_sqr_8_16 ( mt->prod, x );
74 0 : else if( k==16UL ) bignum_ksqr_16_32( mt->prod, x, mt->kt );
75 0 : else if( k==32UL ) bignum_ksqr_32_64( mt->prod, x, mt->kt );
76 0 : else bignum_sqr( 2UL*k, mt->prod, k, x );
77 0 : ulong cy = bignum_emontredc_8n( k, mt->prod, mt->m, mt->w );
78 0 : bignum_optsub( k, out, mt->prod+k, cy, mt->m );
79 0 : return;
80 0 : }
81 0 : #endif
82 0 : bignum_amontsqr( k, out, x, mt->m );
83 0 : }
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 0 : ulong const * y ) {
90 0 : ulong k = mt->k;
91 0 : #if FD_BIG_MOD_EXP_ADX
92 0 : if( FD_LIKELY( mt->fast ) ) {
93 0 : if ( k== 8UL ) bignum_mul_8_16 ( mt->prod, x, y );
94 0 : else if( k==16UL ) bignum_kmul_16_32( mt->prod, x, y, mt->kt );
95 0 : else if( k==32UL ) bignum_kmul_32_64( mt->prod, x, y, mt->kt );
96 0 : else bignum_mul( 2UL*k, mt->prod, k, x, k, y );
97 0 : ulong cy = bignum_emontredc_8n( k, mt->prod, mt->m, mt->w );
98 0 : bignum_optsub( k, out, mt->prod+k, cy, mt->m );
99 0 : return;
100 0 : }
101 0 : #endif
102 0 : bignum_amontmul( k, out, x, y, mt->m );
103 0 : }
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 0 : ulong k ) {
112 0 : ulong km = k; while( !m[km-1UL] ) km--; /* m>1, so km>=1 */
113 0 : ulong ka = k; while( ka && !a[ka-1UL] ) ka--;
114 :
115 0 : fd_memset( r, 0, k*sizeof(ulong) );
116 :
117 : /* a < m: nothing to reduce */
118 0 : if( ka<km ) { fd_memcpy( r, a, ka*sizeof(ulong) ); return; }
119 0 : if( ka==km ) {
120 0 : ulong i = km;
121 0 : while( i && a[i-1UL]==m[i-1UL] ) i--;
122 0 : if( i && a[i-1UL]<m[i-1UL] ) { fd_memcpy( r, a, ka*sizeof(ulong) ); return; }
123 0 : }
124 :
125 0 : ulong mn[ FD_BIG_MOD_EXP_MAX_WORDS ]; /* m << s, km words */
126 0 : ulong an[ FD_BIG_MOD_EXP_MAX_WORDS+1UL ]; /* a << s, ka+1 words */
127 :
128 0 : ulong s = 63UL - (ulong)fd_ulong_find_msb( m[km-1UL] );
129 0 : if( s ) {
130 0 : for( ulong i=km; i-->1UL; ) mn[i] = (m[i]<<s) | (m[i-1UL]>>(64UL-s));
131 0 : mn[0] = m[0]<<s;
132 0 : an[ka] = a[ka-1UL]>>(64UL-s);
133 0 : for( ulong i=ka; i-->1UL; ) an[i] = (a[i]<<s) | (a[i-1UL]>>(64UL-s));
134 0 : an[0] = a[0]<<s;
135 0 : } else {
136 0 : fd_memcpy( mn, m, km*sizeof(ulong) );
137 0 : fd_memcpy( an, a, ka*sizeof(ulong) );
138 0 : an[ka] = 0UL;
139 0 : }
140 :
141 0 : if( km==1UL ) {
142 0 : ulong d = mn[0];
143 0 : uint128 rem = 0;
144 0 : for( ulong j=ka+1UL; j-->0UL; ) rem = ( (rem<<64) | an[j] ) % d; /* rem<d, so no overflow */
145 0 : r[0] = (ulong)rem>>s;
146 0 : return;
147 0 : }
148 :
149 0 : ulong dh = mn[km-1UL];
150 0 : ulong dl = mn[km-2UL];
151 0 : for( ulong j=ka-km+1UL; j-->0UL; ) {
152 : /* quotient digit estimate for the window an[j..j+km]; an[j+km]<=dh */
153 0 : ulong qhat;
154 0 : uint128 rhat;
155 0 : if( an[j+km]>=dh ) {
156 0 : qhat = ~0UL;
157 0 : rhat = (uint128)an[j+km-1UL] + dh;
158 0 : } else {
159 0 : uint128 num = ( (uint128)an[j+km]<<64 ) | an[j+km-1UL];
160 0 : qhat = (ulong)( num/dh );
161 0 : rhat = num%dh;
162 0 : }
163 0 : while( !(rhat>>64) && (uint128)qhat*dl > ( (rhat<<64) | an[j+km-2UL] ) ) {
164 0 : qhat--;
165 0 : rhat += dh;
166 0 : }
167 :
168 : /* an[j..j+km] -= qhat*mn */
169 0 : ulong carry = 0UL;
170 0 : ulong borrow = 0UL;
171 0 : for( ulong i=0UL; i<km; i++ ) {
172 0 : uint128 p = (uint128)qhat*mn[i] + carry;
173 0 : carry = (ulong)( p>>64 );
174 0 : uint128 d = (uint128)an[j+i] - (ulong)p - borrow;
175 0 : an[j+i] = (ulong)d;
176 0 : borrow = (ulong)( d>>64 ) & 1UL;
177 0 : }
178 0 : uint128 d = (uint128)an[j+km] - carry - borrow;
179 0 : an[j+km] = (ulong)d;
180 :
181 : /* estimate one too large: add mn back */
182 0 : if( d>>64 ) {
183 0 : ulong c = 0UL;
184 0 : for( ulong i=0UL; i<km; i++ ) {
185 0 : uint128 t = (uint128)an[j+i] + mn[i] + c;
186 0 : an[j+i] = (ulong)t;
187 0 : c = (ulong)( t>>64 );
188 0 : }
189 0 : an[j+km] += c;
190 0 : }
191 0 : }
192 :
193 : /* remainder is an[0..km), with an[km]==0 */
194 0 : if( s ) {
195 0 : for( ulong i=0UL; i<km; i++ ) r[i] = (an[i]>>s) | (an[i+1UL]<<(64UL-s));
196 0 : } else {
197 0 : fd_memcpy( r, an, km*sizeof(ulong) );
198 0 : }
199 0 : }
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 0 : uchar const * mod, ulong mod_len ) {
206 :
207 0 : if( FD_UNLIKELY( !fd_big_mod_exp_modulus_is_valid( mod, mod_len ) ) )
208 0 : return FD_BIG_MOD_EXP_ERR_MODULUS;
209 :
210 0 : ulong k = ( fd_ulong_max( base_len, mod_len ) + 7UL )/8UL; /* [1,MAX_WORDS] */
211 0 : ulong e = ( exp_len + 7UL )/8UL; /* [0,MAX_WORDS] */
212 :
213 0 : ulong a [ FD_BIG_MOD_EXP_MAX_WORDS ]; /* base, zero extended to k */
214 0 : ulong p [ FD_BIG_MOD_EXP_MAX_WORDS ]; /* exponent, zero extended to e */
215 0 : ulong m [ FD_BIG_MOD_EXP_MAX_WORDS ]; /* modulus, zero extended to k */
216 0 : ulong r2 [ FD_BIG_MOD_EXP_MAX_WORDS ]; /* 2^{128k} mod m */
217 0 : ulong u [ FD_BIG_MOD_EXP_MAX_WORDS ]; /* accumulator, ping */
218 0 : ulong v [ FD_BIG_MOD_EXP_MAX_WORDS ]; /* accumulator, pong */
219 0 : ulong tbl[ FD_BIG_MOD_EXP_TBL_CNT ][ FD_BIG_MOD_EXP_MAX_WORDS ]; /* mont(base^i) */
220 0 : ulong prd[ 2UL*FD_BIG_MOD_EXP_MAX_WORDS ]; /* fast path: product */
221 0 : ulong kt [ FD_BIG_MOD_EXP_KT_WORDS ]; /* fast path: k{sqr,mul} temp */
222 :
223 0 : fd_memset( a, 0, k*sizeof(ulong) );
224 0 : fd_memset( m, 0, k*sizeof(ulong) );
225 0 : fd_memcpy( (uchar *)a, base, base_len );
226 0 : fd_memcpy( (uchar *)m, mod, mod_len );
227 0 : if( FD_LIKELY( e ) ) {
228 0 : fd_memset( p, 0, e*sizeof(ulong) );
229 0 : fd_memcpy( (uchar *)p, exp, exp_len );
230 0 : }
231 :
232 0 : 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 0 : if( FD_UNLIKELY( !nbits ) ) {
236 0 : fd_memset( out, 0, mod_len );
237 0 : out[0] = (uchar)1;
238 0 : return FD_BIG_MOD_EXP_SUCCESS;
239 0 : }
240 :
241 : /* base^1 == base mod m: one reduction, no Montgomery setup. */
242 0 : if( FD_UNLIKELY( nbits==1UL ) ) {
243 0 : fd_big_mod_exp_reduce( u, a, m, k );
244 0 : fd_memcpy( out, (uchar const *)u, mod_len );
245 0 : return FD_BIG_MOD_EXP_SUCCESS;
246 0 : }
247 :
248 : /* Fixed window. */
249 0 : ulong w = 1UL;
250 0 : ulong best = nbits + nbits; /* cost at w==1 */
251 0 : for( ulong t=2UL; t<=FD_BIG_MOD_EXP_MAX_WINDOW; t++ ) {
252 0 : ulong cost = ((1UL<<t)-2UL) + nbits + (nbits+t-1UL)/t;
253 0 : if( cost<best ) { best = cost; w = t; }
254 0 : }
255 :
256 0 : fd_big_mod_exp_mont_t mt[1];
257 0 : mt->m = m;
258 0 : mt->k = k;
259 0 : mt->prod = prd;
260 0 : mt->kt = kt;
261 0 : mt->fast = FD_BIG_MOD_EXP_ADX && k>=8UL && !(k & 7UL);
262 0 : mt->w = mt->fast ? word_negmodinv( m[0] ) : 0UL;
263 :
264 0 : bignum_amontifier( k, r2, m, v /* scratch, >=k words */ );
265 0 : fd_big_mod_exp_mont_mul( mt, tbl[1], r2, a );
266 0 : for( ulong i=2UL; i<(1UL<<w); i++ ) {
267 0 : fd_big_mod_exp_mont_mul( mt, tbl[i], tbl[i-1], tbl[1] );
268 0 : }
269 :
270 0 : ulong top = nbits % w;
271 0 : if( !top ) top = w;
272 0 : ulong pos = nbits - top; /* a multiple of w */
273 :
274 0 : ulong * cur = u;
275 0 : ulong * nxt = v;
276 0 : fd_memcpy( cur, tbl[ fd_big_mod_exp_bits( p, pos, top ) ], k*sizeof(ulong) );
277 :
278 0 : while( pos ) {
279 0 : pos -= w;
280 0 : for( ulong j=0UL; j<w; j++ ) {
281 0 : fd_big_mod_exp_mont_sqr( mt, nxt, cur );
282 0 : ulong * s = cur; cur = nxt; nxt = s;
283 0 : }
284 0 : ulong d = fd_big_mod_exp_bits( p, pos, w );
285 0 : if( d ) {
286 0 : fd_big_mod_exp_mont_mul( mt, nxt, cur, tbl[d] );
287 0 : ulong * s = cur; cur = nxt; nxt = s;
288 0 : }
289 0 : }
290 :
291 0 : bignum_demont( k, nxt, cur, m );
292 0 : fd_memcpy( out, (uchar const *)nxt, mod_len );
293 :
294 0 : return FD_BIG_MOD_EXP_SUCCESS;
295 0 : }
296 :
297 : #endif /* FD_HAS_S2NBIGNUM */
|
