Line data Source code
1 : #ifndef HEADER_fd_src_ballet_bigint_fd_uint256_h
2 : #error "Do not include this directly; use fd_uint256.h"
3 : #endif
4 :
5 : /* Implementation of uint256 Montgomery mul.
6 :
7 : TODOs:
8 : - efficient sqr, Alg. 5 https://eprint.iacr.org/2022/1400
9 : - regular CIOS (the current impl *probably* won't work for secp256k1/r1)
10 :
11 : This is used under the hood to implement:
12 : - bn254 field arithmetic
13 : - bn254 scalar field arithmetic
14 : - (TODO) ed25519 scalar field arithmetic
15 : - (TODO) secp256k1 field and scalar arithmetic
16 : - (TODO) secp256r1 field and scalar arithmetic
17 : - ...
18 : */
19 :
20 : #define INLINE static inline __attribute__((always_inline))
21 :
22 : #ifdef FD_USING_GCC
23 : #define OPTIMIZE __attribute__((optimize("unroll-loops")))
24 : #else
25 : #define OPTIMIZE
26 : #endif
27 :
28 : #if FD_HAS_X86
29 : #if defined(__GNUC__) && !defined(__clang__)
30 : #include <x86gprintrin.h>
31 : #else
32 : #include <immintrin.h>
33 : #endif
34 : #endif
35 :
36 : /* Utility functions for fd_uint256_mul_mod_p.
37 : Implementation is based on uint128.
38 : The implementations WITHOUT uint128 are just for completeness, we
39 : should avoid using them (and, e.g., rely on fiat-crypto mul instead). */
40 :
41 : INLINE void
42 : fd_ulong_sub_borrow(
43 : ulong * r, /* out (r=a-b) */
44 : int * b, /* out borrow flag */
45 : ulong a0,
46 : ulong a1,
47 : int bi /* in borrow flag */
48 204628191 : ) {
49 204628191 : # if FD_HAS_X86
50 204628191 : *b = (uchar)_subborrow_u64( (uchar)bi, a0, a1, (unsigned long long *)r );
51 : # else
52 : a1 += !!bi;
53 : *r = a0 - a1;
54 : *b = a0 < a1;
55 : # endif
56 204628191 : }
57 :
58 : INLINE void
59 : fd_ulong_add_carry(
60 : ulong * r, /* out (r=a+b+ci) */
61 : int * c, /* out carry flag */
62 : ulong a0,
63 : ulong a1,
64 : int ci /* in carry flag */
65 215773138 : ) {
66 215773138 : # if FD_HAS_X86
67 215773138 : *c = (uchar)_addcarry_u64( (uchar)ci, a0, a1, (unsigned long long *)r );
68 : # else
69 : ulong r0 = a0 + a1;
70 : ulong r1 = r0 + !!ci;
71 : *r = r1;
72 : *c = (int)((r0 < a0) | (r1 < r0));
73 : # endif
74 215773138 : }
75 :
76 : #if FD_HAS_INT128
77 :
78 : INLINE void
79 566137600 : fd_ulong_vec_mul( ulong l[4], ulong h[4], ulong const a[4], ulong b ) {
80 566137600 : uint128 r0 = ((uint128)a[0]) * ((uint128)b);
81 566137600 : uint128 r1 = ((uint128)a[1]) * ((uint128)b);
82 566137600 : uint128 r2 = ((uint128)a[2]) * ((uint128)b);
83 566137600 : uint128 r3 = ((uint128)a[3]) * ((uint128)b);
84 566137600 : l[0] = (ulong)r0; h[0] = (ulong)(r0 >> 64);
85 566137600 : l[1] = (ulong)r1; h[1] = (ulong)(r1 >> 64);
86 566137600 : l[2] = (ulong)r2; h[2] = (ulong)(r2 >> 64);
87 566137600 : l[3] = (ulong)r3; h[3] = (ulong)(r3 >> 64);
88 566137600 : }
89 :
90 : INLINE void
91 2265472168 : fd_ulong_add_carry4( ulong *l, uchar *h, ulong a0, ulong a1, ulong a2, uchar a3 ) {
92 2265472168 : uint128 r = ((uint128)a0) + ((uint128)a1) + ((uint128)a2) + ((uint128)a3);
93 2265472168 : *l = (ulong)r;
94 2265472168 : *h = (uchar)(r >> 64);
95 2265472168 : }
96 :
97 : #else
98 :
99 : INLINE void
100 : fd_ulong_mul128( ulong * l, ulong * h, ulong const a, ulong const b ) {
101 : /* First calculate all of the cross products. */
102 : ulong lo_lo = (a & 0xFFFFFFFF) * (b & 0xFFFFFFFF);
103 : ulong hi_lo = (a >> 32) * (b & 0xFFFFFFFF);
104 : ulong lo_hi = (a & 0xFFFFFFFF) * (b >> 32);
105 : ulong hi_hi = (a >> 32) * (b >> 32);
106 :
107 : /* Now add the products together. These will never overflow. */
108 : ulong cross = (lo_lo >> 32) + (hi_lo & 0xFFFFFFFF) + lo_hi;
109 : ulong upper = (hi_lo >> 32) + (cross >> 32) + hi_hi;
110 :
111 : *h = upper;
112 : *l = (cross << 32) | (lo_lo & 0xFFFFFFFF);
113 : }
114 :
115 : INLINE void
116 : fd_ulong_vec_mul( ulong l[4], ulong h[4], ulong const a[4], ulong b ) {
117 : fd_ulong_mul128( &l[0], &h[0], a[0], b );
118 : fd_ulong_mul128( &l[1], &h[1], a[1], b );
119 : fd_ulong_mul128( &l[2], &h[2], a[2], b );
120 : fd_ulong_mul128( &l[3], &h[3], a[3], b );
121 : }
122 :
123 : INLINE void
124 : fd_ulong_add_carry4( ulong *l, uchar *h, ulong a0, ulong a1, ulong a2, uchar a3 ) {
125 : ulong r0 = a0 + a1;
126 : uchar c0 = r0 < a0;
127 :
128 : ulong r1 = a2 + a3;
129 : uchar c1 = r1 < a2;
130 :
131 : *l = r0 + r1;
132 : *h = (uchar)((*l < r0) + c0 + c1);
133 : }
134 :
135 : #endif
136 :
137 : INLINE fd_uint256_t *
138 : fd_uint256_add(fd_uint256_t * r,
139 : fd_uint256_t const * a,
140 230442 : fd_uint256_t const * b ) {
141 230442 : uchar c0;
142 230442 : fd_ulong_add_carry4( &r->limbs[0], &c0, a->limbs[0], b->limbs[0], 0, 0 );
143 230442 : fd_ulong_add_carry4( &r->limbs[1], &c0, a->limbs[1], b->limbs[1], 0, c0 );
144 230442 : fd_ulong_add_carry4( &r->limbs[2], &c0, a->limbs[2], b->limbs[2], 0, c0 );
145 230442 : fd_ulong_add_carry4( &r->limbs[3], &c0, a->limbs[3], b->limbs[3], 0, c0 );
146 230442 : return r;
147 230442 : }
148 :
149 : /* fd_uint256_mul_mod_p computes r = a * b mod p, using the CIOS method.
150 : r, a, b are in Montgomery representation (p is not).
151 :
152 : This is an efficient implementation of CIOS that works when (circa) p < 2^255
153 : (precisely, p->limbs[3] < (2^64-1)/2 - 1).
154 : Alg. 2, https://eprint.iacr.org/2022/1400
155 : Code example, for bn254: https://github.com/Consensys/gnark-crypto/blob/v0.12.1/ecc/bn254/fp/element_ops_purego.go#L66
156 :
157 : In go lang, bits.Add64 has carry 0, 1.
158 : We allow the carry to be a uchar, so we can dp a single add chain after each mul.
159 :
160 : This function is intended to be wrapped into a fd_<field>_mul( r, a, b ).
161 : Experimentally we found that:
162 : 1. We have to force inlining for this function, otherwise compilers tend to reuse
163 : the function, introducing overhead.
164 : 2. In GCC, we have to force loop unrolling optimization *in the outer fd_<field>_mul()*
165 : function, otherwise performance degrades significantly.
166 : For this we added the macro FD_UINT256_FP_MUL_IMPL. */
167 :
168 : #if FD_HAS_X86 && defined(__BMI2__) && defined(__ADX__)
169 : __asm__( ".include \"src/ballet/bigint/fd_uint256_mul.inc\"" );
170 : #endif
171 :
172 : INLINE fd_uint256_t *
173 : fd_uint256_mul_mod_p( fd_uint256_t * r,
174 : fd_uint256_t const * a,
175 : fd_uint256_t const * b,
176 : fd_uint256_t const * p,
177 106150800 : ulong const p_inv ) {
178 35383600 : #if FD_HAS_X86 && defined(__BMI2__) && defined(__ADX__)
179 35383600 : register fd_uint256_t * _r __asm__("rdi") = r;
180 35383600 : register fd_uint256_t const * _a __asm__("rsi") = a;
181 35383600 : register fd_uint256_t const * _b __asm__("rdx") = b;
182 35383600 : register fd_uint256_t const * _p __asm__("rcx") = p;
183 35383600 : register ulong _pinv __asm__("r8") = p_inv;
184 35383600 : __asm__ __volatile__ (
185 35383600 : "_fd_uint256_mul_mod_p %[r], %[a], %[b], %[p], %[pinv]"
186 35383600 : : [r]"+r"(_r), [a]"+r"(_a), [b]"+r"(_b), [pinv]"+r"(_pinv)
187 35383600 : : [p]"r"(_p)
188 35383600 : : "rax", "r9", "r10", "r11", "cc", "memory"
189 35383600 : );
190 35383600 : r = _r;
191 : #else
192 70767200 : ulong FD_ALIGNED t[4] = { 0 };
193 70767200 : ulong FD_ALIGNED u[4];
194 70767200 : ulong FD_ALIGNED h[4];
195 70767200 : ulong FD_ALIGNED l[4];
196 70767200 : uchar c0, c1;
197 70767200 : ulong tmp;
198 70767200 : ulong m;
199 :
200 353836000 : for( int i=0; i<4; i++ ) {
201 283068800 : fd_ulong_vec_mul( l, u, a->limbs, b->limbs[i] );
202 283068800 : fd_ulong_add_carry4( &t[0], &c0, t[0], l[0], 0, 0 );
203 283068800 : fd_ulong_add_carry4( &t[1], &c0, t[1], l[1], u[0], c0 );
204 283068800 : fd_ulong_add_carry4( &t[2], &c0, t[2], l[2], u[1], c0 );
205 283068800 : fd_ulong_add_carry4( &t[3], &c0, t[3], l[3], u[2], c0 );
206 :
207 283068800 : m = t[0] * p_inv;
208 :
209 283068800 : fd_ulong_vec_mul( l, h, p->limbs, m );
210 283068800 : fd_ulong_add_carry4( &tmp, &c1, t[0], l[0], 0, 0 );
211 283068800 : fd_ulong_add_carry4( &t[0], &c1, t[1], l[1], h[0], c1 );
212 283068800 : fd_ulong_add_carry4( &t[1], &c1, t[2], l[2], h[1], c1 );
213 283068800 : fd_ulong_add_carry4( &t[2], &c1, t[3], l[3], h[2], c1 );
214 283068800 : t[3] = u[3] + h[3] + c0 + c1;
215 283068800 : }
216 :
217 70767200 : r->limbs[0] = t[0];
218 70767200 : r->limbs[1] = t[1];
219 70767200 : r->limbs[2] = t[2];
220 70767200 : r->limbs[3] = t[3];
221 :
222 70767200 : if( fd_uint256_cmp( r, p ) >= 0 ) {
223 3620114 : int b = 0;
224 3620114 : fd_ulong_sub_borrow( &r->limbs[0], &b, r->limbs[0], p->limbs[0], b );
225 3620114 : fd_ulong_sub_borrow( &r->limbs[1], &b, r->limbs[1], p->limbs[1], b );
226 3620114 : fd_ulong_sub_borrow( &r->limbs[2], &b, r->limbs[2], p->limbs[2], b );
227 3620114 : fd_ulong_sub_borrow( &r->limbs[3], &b, r->limbs[3], p->limbs[3], b );
228 3620114 : }
229 70767200 : #endif
230 106150800 : return r;
231 106150800 : }
232 :
233 : /* FD_UINT256_FP_MUL_IMPL macro to properly implement Fp mul based on
234 : fd_uint256_mul_mod_p().
235 : In GCC we need to explicitly force loop unroll. */
236 : #define FD_UINT256_FP_MUL_IMPL(fp, p, p_inv) \
237 : static inline fp ## _t * OPTIMIZE \
238 : fp ## _mul( fp ## _t * r, \
239 : fp ## _t const * a, \
240 106150800 : fp ## _t const * b ) { \
241 106150800 : return fd_uint256_mul_mod_p( r, a, b, p, p_inv ); \
242 106150800 : }
|
