• Reasercher and Software Engineer at Cybozu Labs, Inc. • Developer of some OSS • mcl/bls (https://github.com/herumi/mcl) • Pairing/BLS Signature Library • For DFINITY (2016), for Ethereum 2.0 (2019) • https://www.npmjs.com/package/mcl-wasm WeAssembly package.130k DLs/week • Xbyak/Xbyak_aarch64 • https://github.com/herumi/xbyak • JIT assembler for x64/AArch64 (A64) by C++ • Designed For Low-Level Optimization • Heavily Used in oneDNN (Intel AI Framework) / ZenDNN (AMD), Fugaku (Japanese Supercomputer) My Introduction 2 /32 @herumi
• Motivation • For security reasons, we want to write in a "safe language" as much as possible. • But high performance is also required. • We want to minimize the use of assembly language. • What are some safe ways to write such code? • How practical is a generic description using LLVM? • DSL (domain specific language) makes it easier to read and write. • The current performance and limitations of these approaches. Today's Topics 3 /32
• Pairing is a mathematical function used in BLS signature, zk-SNARK, homomorphic encryption, etc. • Multiplication on finite fields occupies most of the time for pairing calculations. • 256~384-bit prime of a finite field is often used. Pairing 4 /32
• Compare with blst (https://github.com/supranational/blst) • Finite field Operations on MacBook Pro, M1, 2020 nsec • DSL is almost at blst speed. • https://github.com/herumi/mcl-ff (under construction) • See below for details. Benchmark of Generated Code by DSL 381 bit prime blst DSL+LLVM add 6.40 6.50 sub 5.30 5.60 mont (mul) 31.40 31.50 5 /32
• A sum of 64-bit registers is 65 bits. • 65bit = 64bit + CF (1bit) • How CF is calculated Review of Addition of Two Digits 36 + 47 3 6 + 4 7 ----- 1 3 3 4 ----- 8 3 increase in digits bool hasCF(U x, U y) { U z = x + y; return z < x; // CF = 1 if z < x } 7 /32
• RISC-V has no carry operations, so this is good, but slow on x64 (x86-64) / A64 (AArch64). addT using C++ // U = uint64_t or uint32_t // ⟨𝑧⟩𝑁 = ⟨𝑥⟩𝑁 + ⟨𝑦⟩𝑁 and return CF template U addT(U *z, const U *x, const U *y) { U c = 0; // There is no CF at first. for (size_t i = 0; i < N; i++) { U xc = x[i] + c; c = xc < c; // CF1 U yi = y[i]; xc += yi; c += xc < yi; // CF2 z[i] = xc; } return c; } 8 /32
• ll (LLVM IR) seems good, but it is hard to write by hand. • I developed a simple DSL to make it easier to write. • Variables can be rewritten. • non SSA • No need to write register size • ll : %x1 = add i256 %x, %y DSL : x = add(x, y); DSL for LLVM using Python (and C++) code (ll) clang –S –O2 x64 A64 RISC-V MIPS generator DSL WASM 16 /32
• gen_add(N) generates the previous code for each N Generating add.ll with DSL unit = 64 def gen_add(N): bit = unit * N pz = IntPtr(unit) # define pointer to int px = IntPtr(unit) py = IntPtr(unit) name = f'mcl_fp_addPre{N}' with Function(name, Int(unit), pz, px, py): x = loadN(px, N) # x = ⟨∗ 𝑝𝑥⟩𝑁 y = loadN(py, N) # y = ⟨∗ 𝑝𝑦⟩𝑁 x = zext(x, bit + unit) y = zext(y, bit + unit) z = add(x, y) storeN(trunc(z, bit), pz) r = trunc(lshr(z, bit), unit) ret(r) 17 /32
• Python code of add of a finite field • Introduce select function Fp::add Using Python # assume 0 <= x, y < p def fp_add(x, y): z = x + y if z >= p: z -= p return z def select(cond, x, y): if cond: return x else: return y def fp_add(x, y): z = x + y w = z – p return select(w < 0, z, w) 18 /32
• gen_fp_add(N) generates F::add for each N Fp::add Using DSL def gen_fp_add(name, N, pp): bit = unit * N pz = IntPtr(unit) px = IntPtr(unit) py = IntPtr(unit) with Function(name, Void, pz, px, py): x = zext(loadN(px, N), bit + unit) y = zext(loadN(py, N), bit + unit) x = add(x, y) # x ← x + y p = zext(load(pp), bit + unit) y = sub(x, p) # y ← x + y - p c = trunc(lshr(y, bit), 1) # c = MSB of y x = select(c, x, y) # x = (y < 0) ? x : y x = trunc(x, bit) storeN(x, pz) Differences with gen_add() 19 /32
• Full-Bit Prime • The bit size of p is a multiple of 64. e.g. 𝑝 = 256, 384 • Non-Full-Bit Prime • e.g. 𝑝 = 381 • If p is non-full-bit, the x+y value does not have to be zext. Optimizing Non-Full-Bit Primes x = loadN(px, N) y = loadN(py, N) x = add(x, y) p = load(pp) y = sub(x, p) c = trunc(lshr(y, bit - 1), 1) # c = MSB(y) = (y < 0) x = select(c, x, y) storeN(x, pz) ret(Void) 20 /32
• u64 mulUnit(u64 *z, const u64 *x, u64 y); • computes 𝑥 𝑁 × 𝑦 = 𝑧 𝑁+1 . • LLVM IR has only u128 mul(u64, u64); • pack is a combination of shl+or Multiplication by u64 [x3:x2:x1:x0] X y ---------------- [H0:L0] [H1:L1] [H2:L2] [H3:L3] ----------------- [z4:z3:z2:z1:z0] L=[ 0:L3:L2:L1:L0] ＋ H=[H3:H2:H1:H0: 0] ↓ z=[z4:z3:z2:z1:z0] with Function(f'mulUnit{N}', z, px, y) Ls = [] Hs = [] y = zext(y, unit*2) for i in range(N): x = load(getelementptr(px, i)) # x[i] xy = mul(zext(x, unit*2), y) Ls.append(trunc(xy, unit)) Hs.append(trunc(lshr(xy, unit), unit)) bu = bit + unit L = zext(pack(Ls), bu) H = shl(zext(pack(Hs), bu), unit) z = add(L, H) 22 /32
• A method for fast modular multiplication • toM(x) : converts x to a Montgomery form • fromM(x) : converts x to a normal form • 𝑥 × 𝑦 ≡ 𝑓𝑟𝑜𝑚𝑀 𝑚𝑜𝑛𝑡 𝑡𝑜𝑀 𝑥 , 𝑡𝑜𝑀 𝑦 mod p • Preparation • 𝐿 = 64 𝑜𝑟 32 and 𝑀 = 2𝐿 • 𝑝 = 𝑝 𝑁 : a prime • 𝑖𝑝 ≡ 𝑝−1 mod 𝑀 Montgomery Multiplication 26 /32
• An example code of mont(x, y) Mont using Python L=64 M=2**L # ip = p^(-1) mod M def mont(x, y): t = 0 for i in range(N): t += x * ((y >> (L * i)) % M) # mulUnit q = ((t % M) * ip) % M # u64 x u64 → u64 t += q * p # mulUnit t >>= L if t >= p: t -= p return t 27 /32
• Assume Non-Full-Bit prime Mont using DSL def gen_mont(name, mont, mulUnit): pz = IntPtr(unit) px = IntPtr(unit) py = IntPtr(unit) with Function(name, Void, pz, px, py): t = Imm(0, unit*N+unit) for i in range(N): y = load(getelementptr(py, i)) t = add(t, call(mulUnit, px, y)) q = mul(trunc(t, unit), ip) t = add(t, call(mulUnit, pp, q)) t = lshr(t, unit) t = trunc(t, unit*N) vc = sub(t, loadN(pp, N)) c = trunc(lshr(vc, unit*N-1), 1) storeN(select(c, t, vc), pz) ret(Void) # Python code t += x * ((y >> (L*i)) % M) q = ((t % M) * ip) % M t += q * p t >>= L if t >= p: t -= p 𝑧 𝑁+1 =mulUnit( 𝑥 𝑁 ,y) 28 /32
• Compare with blst (https://github.com/supranational/blst) • on MacBook Pro, M1, 2020 nsec • mont on Xeon Platinum 8280 (turbo boost off) • DSL+LLVM with –O2 is not fast. • slower than mcl even if –mbmi2 (using mulx) option is added. • mcl (and blst) uses adox/adcx that can be used after IvyBridge. Benchmark of Generated Code by DSL 381 bit prime blst DSL+LLVM add 6.40 6.50 sub 5.30 5.60 mont 31.40 31.50 381 bit prime mcl DSL+LLVM DSL+LLVM –mbmi2 mont 28.93 42.86 35.52 29 /32
• We only have one CF, then z and xy must be kept in registers until they are added. • A64 has 29 regs. but x64 has only 15 regs. • adox/adcx offer two CFs. Register spills can be avoided. • clang/gcc will not use the instruction set now. • https://gcc.gnu.org/legacy-ml/gcc-help/2017-08/msg00100.html Computing 𝑧 𝑁 + 𝑥 𝑁 ∗ 𝑦 with adox/adcx [ z ] [H0 L0] ← 𝑥0 𝑦 [H1 L1] ← 𝑥1 𝑦 [H2 L2] ← 𝑥2 𝑦 [H3 L3] ← 𝑥3 𝑦 ------------------ [z4:z3:z2:z1:z0] adcx adox 30 /32
• How to make mulAdd : z += [px] * rdx • z is an array of registers. • px is a pointer to x. • L, H are temporary registers. • rdx is y. • make Mont using mulAdd DSL for x64 # z[n..0] = z[n-1..0] + px[n-1..0] * rdx def mulAdd(z, px): n = len(z)-1 xor_(z[n], z[n]) for i in range(n): mulx(H, L, ptr(px+i*8)) # [H:L] = px[i] * y adox(z[i], L) # z[i] += L with CF if i == n-1: break adcx(z[i + 1], H) # z[i+1] += H with CF' adox(z[n], H) # z[n] += H with CF adc(z[n], 0) # z[n] += CF' 31 /32
• DSL that generates LLVM IR • Advantages • Easy to read and write • Fast enough on M1 • Disadvantages • Not fully optimized for x64 • May be solved by compiler optimization in the future. • Requires x64-specific DSL now. Conclusion 32 /32
herumi
lnit
uhooi
teamlab
kiyoshiro
maroon1st
jpt_corp_official
ogiwarat
_tanako
tomokusaba
seike460
akarin
akitotsukahara
jcasabona
yeseniaperezcruz
shlominoach
eileencodes
lynnandtonic
paulrobertlloyd
chriscoyier
afnizarnur
shpigford
lara
dougneiner