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