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 30866723 : ) {
49 30866723 : # if FD_HAS_X86
50 30866723 : *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 30866723 : }
57 :
58 : #if FD_HAS_INT128
59 :
60 : INLINE void
61 777239360 : fd_ulong_vec_mul( ulong l[4], ulong h[4], ulong const a[4], ulong b ) {
62 777239360 : uint128 r0 = ((uint128)a[0]) * ((uint128)b);
63 777239360 : uint128 r1 = ((uint128)a[1]) * ((uint128)b);
64 777239360 : uint128 r2 = ((uint128)a[2]) * ((uint128)b);
65 777239360 : uint128 r3 = ((uint128)a[3]) * ((uint128)b);
66 777239360 : l[0] = (ulong)r0; h[0] = (ulong)(r0 >> 64);
67 777239360 : l[1] = (ulong)r1; h[1] = (ulong)(r1 >> 64);
68 777239360 : l[2] = (ulong)r2; h[2] = (ulong)(r2 >> 64);
69 777239360 : l[3] = (ulong)r3; h[3] = (ulong)(r3 >> 64);
70 777239360 : }
71 :
72 : INLINE void
73 3109879208 : fd_ulong_add_carry4( ulong *l, uchar *h, ulong a0, ulong a1, ulong a2, uchar a3 ) {
74 3109879208 : uint128 r = ((uint128)a0) + ((uint128)a1) + ((uint128)a2) + ((uint128)a3);
75 3109879208 : *l = (ulong)r;
76 3109879208 : *h = (uchar)(r >> 64);
77 3109879208 : }
78 :
79 : #else
80 :
81 : INLINE void
82 : fd_ulong_mul128( ulong * l, ulong * h, ulong const a, ulong const b ) {
83 : /* First calculate all of the cross products. */
84 : ulong lo_lo = (a & 0xFFFFFFFF) * (b & 0xFFFFFFFF);
85 : ulong hi_lo = (a >> 32) * (b & 0xFFFFFFFF);
86 : ulong lo_hi = (a & 0xFFFFFFFF) * (b >> 32);
87 : ulong hi_hi = (a >> 32) * (b >> 32);
88 :
89 : /* Now add the products together. These will never overflow. */
90 : ulong cross = (lo_lo >> 32) + (hi_lo & 0xFFFFFFFF) + lo_hi;
91 : ulong upper = (hi_lo >> 32) + (cross >> 32) + hi_hi;
92 :
93 : *h = upper;
94 : *l = (cross << 32) | (lo_lo & 0xFFFFFFFF);
95 : }
96 :
97 : INLINE void
98 : fd_ulong_vec_mul( ulong l[4], ulong h[4], ulong const a[4], ulong b ) {
99 : fd_ulong_mul128( &l[0], &h[0], a[0], b );
100 : fd_ulong_mul128( &l[1], &h[1], a[1], b );
101 : fd_ulong_mul128( &l[2], &h[2], a[2], b );
102 : fd_ulong_mul128( &l[3], &h[3], a[3], b );
103 : }
104 :
105 : INLINE void
106 : fd_ulong_add_carry4( ulong *l, uchar *h, ulong a0, ulong a1, ulong a2, uchar a3 ) {
107 : ulong r0 = a0 + a1;
108 : uchar c0 = r0 < a0;
109 :
110 : ulong r1 = a2 + a3;
111 : uchar c1 = r1 < a2;
112 :
113 : *l = r0 + r1;
114 : *h = (uchar)((*l < r0) + c0 + c1);
115 : }
116 :
117 : #endif
118 :
119 : INLINE fd_uint256_t *
120 : fd_uint256_add(fd_uint256_t * r,
121 : fd_uint256_t const * a,
122 230442 : fd_uint256_t const * b ) {
123 230442 : uchar c0;
124 230442 : fd_ulong_add_carry4( &r->limbs[0], &c0, a->limbs[0], b->limbs[0], 0, 0 );
125 230442 : fd_ulong_add_carry4( &r->limbs[1], &c0, a->limbs[1], b->limbs[1], 0, c0 );
126 230442 : fd_ulong_add_carry4( &r->limbs[2], &c0, a->limbs[2], b->limbs[2], 0, c0 );
127 230442 : fd_ulong_add_carry4( &r->limbs[3], &c0, a->limbs[3], b->limbs[3], 0, c0 );
128 230442 : return r;
129 230442 : }
130 :
131 : /* fd_uint256_mul_mod_p computes r = a * b mod p, using the CIOS method.
132 : r, a, b are in Montgomery representation (p is not).
133 :
134 : This is an efficient implementation of CIOS that works when (circa) p < 2^255
135 : (precisely, p->limbs[3] < (2^64-1)/2 - 1).
136 : Alg. 2, https://eprint.iacr.org/2022/1400
137 : Code example, for bn254: https://github.com/Consensys/gnark-crypto/blob/v0.12.1/ecc/bn254/fp/element_ops_purego.go#L66
138 :
139 : In go lang, bits.Add64 has carry 0, 1.
140 : We allow the carry to be a uchar, so we can dp a single add chain after each mul.
141 :
142 : This function is intended to be wrapped into a fd_<field>_mul( r, a, b ).
143 : Experimentally we found that:
144 : 1. We have to force inlining for this function, otherwise compilers tend to reuse
145 : the function, introducing overhead.
146 : 2. In GCC, we have to force loop unrolling optimization *in the outer fd_<field>_mul()*
147 : function, otherwise performance degrades significantly.
148 : For this we added the macro FD_UINT256_FP_MUL_IMPL. */
149 :
150 : #if FD_HAS_X86 && defined(__BMI2__) && defined(__ADX__)
151 : __asm__( ".include \"src/ballet/bigint/fd_uint256_mul.inc\"" );
152 : #endif
153 :
154 : INLINE fd_uint256_t *
155 : fd_uint256_mul_mod_p( fd_uint256_t * r,
156 : fd_uint256_t const * a,
157 : fd_uint256_t const * b,
158 : fd_uint256_t const * p,
159 145732380 : ulong const p_inv ) {
160 48577460 : #if FD_HAS_X86 && defined(__BMI2__) && defined(__ADX__)
161 48577460 : register fd_uint256_t * _r __asm__("rdi") = r;
162 48577460 : register fd_uint256_t const * _a __asm__("rsi") = a;
163 48577460 : register fd_uint256_t const * _b __asm__("rdx") = b;
164 48577460 : register fd_uint256_t const * _p __asm__("rcx") = p;
165 48577460 : register ulong _pinv __asm__("r8") = p_inv;
166 48577460 : __asm__ __volatile__ (
167 48577460 : "_fd_uint256_mul_mod_p %[r], %[a], %[b], %[p], %[pinv]"
168 48577460 : : [r]"+r"(_r), [a]"+r"(_a), [b]"+r"(_b), [pinv]"+r"(_pinv)
169 48577460 : : [p]"r"(_p)
170 48577460 : : "rax", "r9", "r10", "r11", "cc", "memory"
171 48577460 : );
172 48577460 : r = _r;
173 : #else
174 97154920 : ulong FD_ALIGNED t[4] = { 0 };
175 97154920 : ulong FD_ALIGNED u[4];
176 97154920 : ulong FD_ALIGNED h[4];
177 97154920 : ulong FD_ALIGNED l[4];
178 97154920 : uchar c0, c1;
179 97154920 : ulong tmp;
180 97154920 : ulong m;
181 :
182 485774600 : for( int i=0; i<4; i++ ) {
183 388619680 : fd_ulong_vec_mul( l, u, a->limbs, b->limbs[i] );
184 388619680 : fd_ulong_add_carry4( &t[0], &c0, t[0], l[0], 0, 0 );
185 388619680 : fd_ulong_add_carry4( &t[1], &c0, t[1], l[1], u[0], c0 );
186 388619680 : fd_ulong_add_carry4( &t[2], &c0, t[2], l[2], u[1], c0 );
187 388619680 : fd_ulong_add_carry4( &t[3], &c0, t[3], l[3], u[2], c0 );
188 :
189 388619680 : m = t[0] * p_inv;
190 :
191 388619680 : fd_ulong_vec_mul( l, h, p->limbs, m );
192 388619680 : fd_ulong_add_carry4( &tmp, &c1, t[0], l[0], 0, 0 );
193 388619680 : fd_ulong_add_carry4( &t[0], &c1, t[1], l[1], h[0], c1 );
194 388619680 : fd_ulong_add_carry4( &t[1], &c1, t[2], l[2], h[1], c1 );
195 388619680 : fd_ulong_add_carry4( &t[2], &c1, t[3], l[3], h[2], c1 );
196 388619680 : t[3] = u[3] + h[3] + c0 + c1;
197 388619680 : }
198 :
199 97154920 : r->limbs[0] = t[0];
200 97154920 : r->limbs[1] = t[1];
201 97154920 : r->limbs[2] = t[2];
202 97154920 : r->limbs[3] = t[3];
203 :
204 97154920 : if( fd_uint256_cmp( r, p ) >= 0 ) {
205 4716638 : int b = 0;
206 4716638 : fd_ulong_sub_borrow( &r->limbs[0], &b, r->limbs[0], p->limbs[0], b );
207 4716638 : fd_ulong_sub_borrow( &r->limbs[1], &b, r->limbs[1], p->limbs[1], b );
208 4716638 : fd_ulong_sub_borrow( &r->limbs[2], &b, r->limbs[2], p->limbs[2], b );
209 4716638 : fd_ulong_sub_borrow( &r->limbs[3], &b, r->limbs[3], p->limbs[3], b );
210 4716638 : }
211 97154920 : #endif
212 145732380 : return r;
213 145732380 : }
214 :
215 : /* FD_UINT256_FP_MUL_IMPL macro to properly implement Fp mul based on
216 : fd_uint256_mul_mod_p().
217 : In GCC we need to explicitly force loop unroll. */
218 : #define FD_UINT256_FP_MUL_IMPL(fp, p, p_inv) \
219 : static inline fp ## _t * OPTIMIZE \
220 : fp ## _mul( fp ## _t * r, \
221 : fp ## _t const * a, \
222 145732380 : fp ## _t const * b ) { \
223 145732380 : return fd_uint256_mul_mod_p( r, a, b, p, p_inv ); \
224 145732380 : }
|
