commit dc6ccae66ec7a9ef523d4963cdf62a08cb1b50db
parent 541e883fc36426224a149b058fa20b70e195aa13
Author: Jared Tobin <jared@jtobin.io>
Date: Sat, 31 Oct 2026 14:49:27 +0400
Merge branch 'impl/x86'
Adds native Poly1305 on x86_64: the scalar C kernel everywhere, and an
AVX2 4-way kernel where the CPU supports it, alongside the existing
ARM NEON path.
* The internal Crypto.MAC.Poly1305.Arm module and cbits/poly1305_arm.c
are renamed to Crypto.MAC.Poly1305.Native and
cbits/poly1305_native.c, exporting poly1305_mac_native and
poly1305_native_available.
* The scalar 26-bit-limb kernel is portable C, and now also serves
x86_64.
* On x86_64, an AVX2 counterpart of the NEON 4-way kernel handles
64-byte groups from 128 bytes up. AVX2 support (with OS-enabled YMM
state) is checked via CPUID and cached; CPUs without it keep the
scalar kernel.
* A new test compares the dispatching 'mac' against the pure Haskell
loop for every message length from 0 to 640 bytes.
Tested on aarch64-darwin (NEON path) and on an x86_64-linux AMD EPYC
7642 host (AVX2 path), with the LLVM backend, and under ASan/UBSan on
the latter. Under Intel SDE the tests pass emulating Haswell and
Sandy Bridge, the latter falling back to the scalar kernel. On the
x86_64 host, versus the pure Haskell implementation:
mac (114B msg): 347 ns -> 234 ns
mac (1024B msg): 2.33 us -> 572 ns
mac (4096B msg): 8.35 us -> 1.77 us
Diffstat:
8 files changed, 654 insertions(+), 428 deletions(-)
diff --git a/README.md b/README.md
@@ -31,10 +31,11 @@ Haddocks (API documentation, etc.) are hosted at
## Performance
-The aim is best-in-class performance. Current benchmark figures on the
-simple "sunscreen input" from RFC8439 on an M4 Silicon MacBook Air,
-where we avail of hardware acceleration via ARM NEON intrinsics, look
-like (use `cabal bench` to run the benchmark suite):
+The aim is best-in-class performance. We avail of hardware
+acceleration via ARM NEON intrinsics on aarch64, and x86 AVX2
+intrinsics on x86_64 whenever the CPU supports them. Current benchmark
+figures on the simple "sunscreen input" from RFC8439 on an M4 Silicon
+MacBook Air look like (use `cabal bench` to run the benchmark suite):
```
benchmarking ppad-poly1305/mac (big key)
@@ -44,7 +45,7 @@ like (use `cabal bench` to run the benchmark suite):
std dev 742.4 ps (489.7 ps .. 1.169 ns)
```
-On longer inputs the NEON 4-way parallel kernel kicks in, with
+On longer inputs the 4-way parallel kernel kicks in, with
correspondingly better throughput:
```
diff --git a/cbits/poly1305_arm.c b/cbits/poly1305_arm.c
@@ -1,360 +0,0 @@
-#include <stddef.h>
-#include <stdint.h>
-#include <string.h>
-
-#if defined(__aarch64__)
-
-#include <arm_neon.h>
-
-/*
- * Poly1305 (RFC 8439). ARM acceleration via two paths:
- *
- * 1. A scalar 26-bit-limb kernel using 64-bit native arithmetic for
- * the limb multiplications. Used for setup (precomputing r^k),
- * for messages shorter than 64 bytes, and for the < 4-block tail
- * of longer messages. Also serves as the stage-1 reference; the
- * NEON kernel below was added in a separate commit.
- *
- * 2. A NEON 4-way parallel kernel (when at least 4 full blocks
- * remain). Layout: 5 'uint32x4_t' limb vectors hold one limb
- * position each across 4 message blocks; matching r-power
- * vectors hold (r^4, r^3, r^2, r^1) at the same limb position.
- * For each output position d_j we accumulate 5 partial products
- * across both vector halves with 'vmull'/'vmlal', then
- * horizontally sum the 4 lanes with 'vaddvq_u64'. The same
- * carry-propagation pattern as the scalar path reduces back to
- * 26-bit limbs.
- *
- * Reduction rule (used in both paths): 2^130 = 5 (mod 2^130 - 5),
- * so any partial product whose 'position' exceeds 4 folds back as
- * (5 * value) into position (pos - 5).
- */
-
-#define MASK26 0x3ffffffu
-
-/*
- * Multiply two 130-bit values mod (2^130 - 5). Inputs in 5x 26-bit
- * limb form, output in 5x 26-bit limb form (each output limb < 2^26
- * except possibly limb 1, which may carry a small excess absorbed
- * by the next 'mul_mod_p' or by 'normalize').
- */
-static void mul_mod_p(const uint32_t a[5], const uint32_t b[5],
- uint32_t out[5]) {
- uint64_t d0, d1, d2, d3, d4, c;
-
- d0 = (uint64_t)a[0]*b[0]
- + 5 * ((uint64_t)a[4]*b[1] + (uint64_t)a[3]*b[2]
- + (uint64_t)a[2]*b[3] + (uint64_t)a[1]*b[4]);
- d1 = (uint64_t)a[0]*b[1] + (uint64_t)a[1]*b[0]
- + 5 * ((uint64_t)a[4]*b[2] + (uint64_t)a[3]*b[3]
- + (uint64_t)a[2]*b[4]);
- d2 = (uint64_t)a[0]*b[2] + (uint64_t)a[1]*b[1]
- + (uint64_t)a[2]*b[0]
- + 5 * ((uint64_t)a[4]*b[3] + (uint64_t)a[3]*b[4]);
- d3 = (uint64_t)a[0]*b[3] + (uint64_t)a[1]*b[2]
- + (uint64_t)a[2]*b[1] + (uint64_t)a[3]*b[0]
- + 5 * ((uint64_t)a[4]*b[4]);
- d4 = (uint64_t)a[0]*b[4] + (uint64_t)a[1]*b[3]
- + (uint64_t)a[2]*b[2] + (uint64_t)a[3]*b[1]
- + (uint64_t)a[4]*b[0];
-
- c = d0 >> 26; d0 &= MASK26; d1 += c;
- c = d1 >> 26; d1 &= MASK26; d2 += c;
- c = d2 >> 26; d2 &= MASK26; d3 += c;
- c = d3 >> 26; d3 &= MASK26; d4 += c;
- c = d4 >> 26; d4 &= MASK26; d0 += c * 5;
- c = d0 >> 26; d0 &= MASK26; d1 += c;
-
- out[0] = (uint32_t)d0;
- out[1] = (uint32_t)d1;
- out[2] = (uint32_t)d2;
- out[3] = (uint32_t)d3;
- out[4] = (uint32_t)d4;
-}
-
-/*
- * Parse 16 little-endian bytes plus a 'hibit' value (0 or 1) at bit
- * 128 into 5 26-bit limbs.
- */
-static inline void blk2limbs(const uint8_t m[16], uint32_t hibit,
- uint32_t l[5]) {
- uint32_t t0, t1, t2, t3;
- memcpy(&t0, m, 4);
- memcpy(&t1, m + 4, 4);
- memcpy(&t2, m + 8, 4);
- memcpy(&t3, m + 12, 4);
- l[0] = t0 & MASK26;
- l[1] = ((t0 >> 26) | (t1 << 6)) & MASK26;
- l[2] = ((t1 >> 20) | (t2 << 12)) & MASK26;
- l[3] = ((t2 >> 14) | (t3 << 18)) & MASK26;
- l[4] = (t3 >> 8) | (hibit << 24);
-}
-
-/*
- * Process one full 16-byte block: h := (h + m) * r mod p, scalar
- * implementation.
- */
-static inline void scalar_block(uint32_t h[5], const uint32_t r[5],
- const uint8_t m[16], uint32_t hibit) {
- uint32_t blk[5];
- blk2limbs(m, hibit, blk);
-
- uint32_t hl[5];
- hl[0] = h[0] + blk[0];
- hl[1] = h[1] + blk[1];
- hl[2] = h[2] + blk[2];
- hl[3] = h[3] + blk[3];
- hl[4] = h[4] + blk[4];
-
- mul_mod_p(hl, r, h);
-}
-
-/*
- * 4-block NEON update. Computes:
- * h := (h + m_0)*r^4 + m_1*r^3 + m_2*r^2 + m_3*r^1 mod p
- *
- * where m_0..m_3 are four consecutive 16-byte input blocks (so the
- * function consumes 64 bytes). The polynomial identity is the
- * standard Horner expansion of 4 sequential block updates.
- */
-static void neon4_block(uint32_t h[5],
- const uint32_t r1[5], const uint32_t r2[5],
- const uint32_t r3[5], const uint32_t r4[5],
- const uint8_t m[64]) {
- /* Limbify the four blocks. */
- uint32_t b0[5], b1[5], b2[5], b3[5];
- blk2limbs(m, 1, b0);
- blk2limbs(m + 16, 1, b1);
- blk2limbs(m + 32, 1, b2);
- blk2limbs(m + 48, 1, b3);
-
- /* Fold the running accumulator into the first block. */
- b0[0] += h[0]; b0[1] += h[1]; b0[2] += h[2];
- b0[3] += h[3]; b0[4] += h[4];
-
- /* Pack messages: mv[i] = (b0[i], b1[i], b2[i], b3[i]). */
- uint32x4_t mv0 = { b0[0], b1[0], b2[0], b3[0] };
- uint32x4_t mv1 = { b0[1], b1[1], b2[1], b3[1] };
- uint32x4_t mv2 = { b0[2], b1[2], b2[2], b3[2] };
- uint32x4_t mv3 = { b0[3], b1[3], b2[3], b3[3] };
- uint32x4_t mv4 = { b0[4], b1[4], b2[4], b3[4] };
-
- /* Pack r-powers: rv[i] = (r^4[i], r^3[i], r^2[i], r^1[i]). */
- uint32x4_t rv0 = { r4[0], r3[0], r2[0], r1[0] };
- uint32x4_t rv1 = { r4[1], r3[1], r2[1], r1[1] };
- uint32x4_t rv2 = { r4[2], r3[2], r2[2], r1[2] };
- uint32x4_t rv3 = { r4[3], r3[3], r2[3], r1[3] };
- uint32x4_t rv4 = { r4[4], r3[4], r2[4], r1[4] };
-
- /* 5 * r-powers, for partial products whose position wraps past
- * 4 (mod 2^130 - 5 = 5). */
- uint32x4_t rv1_5 = vaddq_u32(vshlq_n_u32(rv1, 2), rv1);
- uint32x4_t rv2_5 = vaddq_u32(vshlq_n_u32(rv2, 2), rv2);
- uint32x4_t rv3_5 = vaddq_u32(vshlq_n_u32(rv3, 2), rv3);
- uint32x4_t rv4_5 = vaddq_u32(vshlq_n_u32(rv4, 2), rv4);
-
- /*
- * Output limb j = sum_i (mv[i] * appropriate r-power for shift j-i,
- * with 5x multiplier when j - i is negative). Per output we sum
- * across the 4 lanes ('vaddvq_u64') after collecting both vmull
- * halves.
- */
-#define MUL5_LO(m_, r_) vmull_u32(vget_low_u32(m_), vget_low_u32(r_))
-#define MUL5_HI(m_, r_) vmull_high_u32(m_, r_)
-#define MLA_LO(acc, m_, r_) \
- vmlal_u32(acc, vget_low_u32(m_), vget_low_u32(r_))
-#define MLA_HI(acc, m_, r_) vmlal_high_u32(acc, m_, r_)
-
- uint64x2_t l0 = MUL5_LO(mv0, rv0);
- uint64x2_t h0 = MUL5_HI(mv0, rv0);
- l0 = MLA_LO(l0, mv1, rv4_5); h0 = MLA_HI(h0, mv1, rv4_5);
- l0 = MLA_LO(l0, mv2, rv3_5); h0 = MLA_HI(h0, mv2, rv3_5);
- l0 = MLA_LO(l0, mv3, rv2_5); h0 = MLA_HI(h0, mv3, rv2_5);
- l0 = MLA_LO(l0, mv4, rv1_5); h0 = MLA_HI(h0, mv4, rv1_5);
- uint64_t d0 = vaddvq_u64(l0) + vaddvq_u64(h0);
-
- uint64x2_t l1 = MUL5_LO(mv0, rv1);
- uint64x2_t h1 = MUL5_HI(mv0, rv1);
- l1 = MLA_LO(l1, mv1, rv0 ); h1 = MLA_HI(h1, mv1, rv0 );
- l1 = MLA_LO(l1, mv2, rv4_5); h1 = MLA_HI(h1, mv2, rv4_5);
- l1 = MLA_LO(l1, mv3, rv3_5); h1 = MLA_HI(h1, mv3, rv3_5);
- l1 = MLA_LO(l1, mv4, rv2_5); h1 = MLA_HI(h1, mv4, rv2_5);
- uint64_t d1 = vaddvq_u64(l1) + vaddvq_u64(h1);
-
- uint64x2_t l2 = MUL5_LO(mv0, rv2);
- uint64x2_t h2 = MUL5_HI(mv0, rv2);
- l2 = MLA_LO(l2, mv1, rv1 ); h2 = MLA_HI(h2, mv1, rv1 );
- l2 = MLA_LO(l2, mv2, rv0 ); h2 = MLA_HI(h2, mv2, rv0 );
- l2 = MLA_LO(l2, mv3, rv4_5); h2 = MLA_HI(h2, mv3, rv4_5);
- l2 = MLA_LO(l2, mv4, rv3_5); h2 = MLA_HI(h2, mv4, rv3_5);
- uint64_t d2 = vaddvq_u64(l2) + vaddvq_u64(h2);
-
- uint64x2_t l3 = MUL5_LO(mv0, rv3);
- uint64x2_t h3 = MUL5_HI(mv0, rv3);
- l3 = MLA_LO(l3, mv1, rv2 ); h3 = MLA_HI(h3, mv1, rv2 );
- l3 = MLA_LO(l3, mv2, rv1 ); h3 = MLA_HI(h3, mv2, rv1 );
- l3 = MLA_LO(l3, mv3, rv0 ); h3 = MLA_HI(h3, mv3, rv0 );
- l3 = MLA_LO(l3, mv4, rv4_5); h3 = MLA_HI(h3, mv4, rv4_5);
- uint64_t d3 = vaddvq_u64(l3) + vaddvq_u64(h3);
-
- uint64x2_t l4 = MUL5_LO(mv0, rv4);
- uint64x2_t h4 = MUL5_HI(mv0, rv4);
- l4 = MLA_LO(l4, mv1, rv3); h4 = MLA_HI(h4, mv1, rv3);
- l4 = MLA_LO(l4, mv2, rv2); h4 = MLA_HI(h4, mv2, rv2);
- l4 = MLA_LO(l4, mv3, rv1); h4 = MLA_HI(h4, mv3, rv1);
- l4 = MLA_LO(l4, mv4, rv0); h4 = MLA_HI(h4, mv4, rv0);
- uint64_t d4 = vaddvq_u64(l4) + vaddvq_u64(h4);
-
-#undef MUL5_LO
-#undef MUL5_HI
-#undef MLA_LO
-#undef MLA_HI
-
- /* Carry propagation, same shape as 'mul_mod_p'. */
- uint64_t c;
- c = d0 >> 26; d0 &= MASK26; d1 += c;
- c = d1 >> 26; d1 &= MASK26; d2 += c;
- c = d2 >> 26; d2 &= MASK26; d3 += c;
- c = d3 >> 26; d3 &= MASK26; d4 += c;
- c = d4 >> 26; d4 &= MASK26; d0 += c * 5;
- c = d0 >> 26; d0 &= MASK26; d1 += c;
-
- h[0] = (uint32_t)d0;
- h[1] = (uint32_t)d1;
- h[2] = (uint32_t)d2;
- h[3] = (uint32_t)d3;
- h[4] = (uint32_t)d4;
-}
-
-/*
- * Compute a 16-byte Poly1305 MAC over 'msg' (length 'msglen') using
- * the 32-byte 'key'. Writes the tag to 'mac_out'.
- */
-void poly1305_mac_arm(const uint8_t key[32], const uint8_t *msg,
- size_t msglen, uint8_t mac_out[16]) {
- /* clamp r */
- uint32_t t0, t1, t2, t3;
- memcpy(&t0, key, 4);
- memcpy(&t1, key + 4, 4);
- memcpy(&t2, key + 8, 4);
- memcpy(&t3, key + 12, 4);
- t0 &= 0x0fffffffu;
- t1 &= 0x0ffffffcu;
- t2 &= 0x0ffffffcu;
- t3 &= 0x0ffffffcu;
-
- uint32_t r[5];
- r[0] = t0 & MASK26;
- r[1] = ((t0 >> 26) | (t1 << 6)) & MASK26;
- r[2] = ((t1 >> 20) | (t2 << 12)) & MASK26;
- r[3] = ((t2 >> 14) | (t3 << 18)) & MASK26;
- r[4] = (t3 >> 8);
-
- uint32_t h[5] = { 0, 0, 0, 0, 0 };
-
- size_t pos = 0;
-
- /* NEON 4-way path: amortizing the r^2/r^3/r^4 precomputation
- * (3 scalar mul_mod_p calls) needs about 4 NEON iterations to
- * break even versus the scalar block loop, so only engage it
- * when we have at least 16 full blocks (256 bytes). */
- if (msglen >= 256) {
- uint32_t r2[5], r3[5], r4[5];
- mul_mod_p(r, r, r2);
- mul_mod_p(r2, r, r3);
- mul_mod_p(r3, r, r4);
-
- while (pos + 64 <= msglen) {
- neon4_block(h, r, r2, r3, r4, msg + pos);
- pos += 64;
- }
- }
-
- /* Scalar tail: any remaining full blocks (< 4 of them). */
- while (pos + 16 <= msglen) {
- scalar_block(h, r, msg + pos, 1);
- pos += 16;
- }
-
- /* Final partial block (1..15 trailing bytes), if any. */
- if (pos < msglen) {
- size_t rem = msglen - pos;
- uint8_t pad[16] = { 0 };
- memcpy(pad, msg + pos, rem);
- pad[rem] = 1;
- scalar_block(h, r, pad, 0);
- }
-
- /* normalize h (mul_mod_p may leave a small excess in h[1]) */
- {
- uint32_t c;
- c = h[1] >> 26; h[1] &= MASK26; h[2] += c;
- c = h[2] >> 26; h[2] &= MASK26; h[3] += c;
- c = h[3] >> 26; h[3] &= MASK26; h[4] += c;
- c = h[4] >> 26; h[4] &= MASK26; h[0] += c * 5;
- c = h[0] >> 26; h[0] &= MASK26; h[1] += c;
- }
-
- /* full reduction to [0, p) via constant-time conditional
- * subtraction of p */
- uint32_t g[5];
- uint32_t c = 5;
- g[0] = h[0] + c; c = g[0] >> 26; g[0] &= MASK26;
- g[1] = h[1] + c; c = g[1] >> 26; g[1] &= MASK26;
- g[2] = h[2] + c; c = g[2] >> 26; g[2] &= MASK26;
- g[3] = h[3] + c; c = g[3] >> 26; g[3] &= MASK26;
- g[4] = h[4] + c;
- uint32_t carry = g[4] >> 26;
- g[4] &= MASK26;
- uint32_t mask = (uint32_t)0 - carry;
- h[0] = (h[0] & ~mask) | (g[0] & mask);
- h[1] = (h[1] & ~mask) | (g[1] & mask);
- h[2] = (h[2] & ~mask) | (g[2] & mask);
- h[3] = (h[3] & ~mask) | (g[3] & mask);
- h[4] = (h[4] & ~mask) | (g[4] & mask);
-
- /* repack 5x 26-bit limbs into 4x 32-bit limbs */
- uint32_t h0 = h[0] | (h[1] << 26);
- uint32_t h1 = (h[1] >> 6) | (h[2] << 20);
- uint32_t h2 = (h[2] >> 12) | (h[3] << 14);
- uint32_t h3 = (h[3] >> 18) | (h[4] << 8);
-
- /* add s (high 16 bytes of key), mod 2^128 */
- uint32_t s0, s1, s2, s3;
- memcpy(&s0, key + 16, 4);
- memcpy(&s1, key + 20, 4);
- memcpy(&s2, key + 24, 4);
- memcpy(&s3, key + 28, 4);
-
- uint64_t a0 = (uint64_t)h0 + s0;
- uint64_t a1 = (uint64_t)h1 + s1 + (a0 >> 32);
- uint64_t a2 = (uint64_t)h2 + s2 + (a1 >> 32);
- uint64_t a3 = (uint64_t)h3 + s3 + (a2 >> 32);
-
- uint32_t o0 = (uint32_t)a0;
- uint32_t o1 = (uint32_t)a1;
- uint32_t o2 = (uint32_t)a2;
- uint32_t o3 = (uint32_t)a3;
-
- memcpy(mac_out + 0, &o0, 4);
- memcpy(mac_out + 4, &o1, 4);
- memcpy(mac_out + 8, &o2, 4);
- memcpy(mac_out + 12, &o3, 4);
-}
-
-int poly1305_arm_available(void) {
- return 1;
-}
-
-#else
-
-void poly1305_mac_arm(const uint8_t *key, const uint8_t *msg,
- size_t msglen, uint8_t *mac_out) {
- (void)key; (void)msg; (void)msglen; (void)mac_out;
-}
-
-int poly1305_arm_available(void) {
- return 0;
-}
-
-#endif
diff --git a/cbits/poly1305_native.c b/cbits/poly1305_native.c
@@ -0,0 +1,554 @@
+#include <stddef.h>
+#include <stdint.h>
+#include <string.h>
+
+#if defined(__aarch64__)
+#include <arm_neon.h>
+#elif defined(__x86_64__)
+#include <cpuid.h>
+#include <immintrin.h>
+#endif
+
+#if defined(__aarch64__) || defined(__x86_64__)
+
+/*
+ * Poly1305 (RFC 8439), via two paths:
+ *
+ * 1. A portable scalar 26-bit-limb kernel using 64-bit native
+ * arithmetic for the limb multiplications, compiled on both
+ * aarch64 and x86_64. Used for setup (precomputing r^k),
+ * for messages shorter than 64 bytes, and for the < 4-block tail
+ * of longer messages. Also serves as the stage-1 reference; the
+ * NEON kernel below was added in a separate commit.
+ *
+ * 2. A 4-way parallel kernel (when at least 4 full blocks remain).
+ * On aarch64 it uses NEON: 5 'uint32x4_t' limb vectors hold one
+ * limb position each across 4 message blocks; matching r-power
+ * vectors hold (r^4, r^3, r^2, r^1) at the same limb position.
+ * For each output position d_j we accumulate 5 partial products
+ * across both vector halves with 'vmull'/'vmlal', then
+ * horizontally sum the 4 lanes with 'vaddvq_u64'. The same
+ * carry-propagation pattern as the scalar path reduces back to
+ * 26-bit limbs.
+ *
+ * On x86_64 it uses AVX2, where the CPU supports it, with the
+ * same structure: one '__m256i' per limb position, one 64-bit
+ * lane per block, and '_mm256_mul_epu32' in place of
+ * 'vmull'/'vmlal'.
+ *
+ * Reduction rule (used in both paths): 2^130 = 5 (mod 2^130 - 5),
+ * so any partial product whose 'position' exceeds 4 folds back as
+ * (5 * value) into position (pos - 5).
+ */
+
+#define MASK26 0x3ffffffu
+
+/*
+ * Multiply two 130-bit values mod (2^130 - 5). Inputs in 5x 26-bit
+ * limb form, output in 5x 26-bit limb form (each output limb < 2^26
+ * except possibly limb 1, which may carry a small excess absorbed
+ * by the next 'mul_mod_p' or by 'normalize').
+ */
+static void mul_mod_p(const uint32_t a[5], const uint32_t b[5],
+ uint32_t out[5]) {
+ uint64_t d0, d1, d2, d3, d4, c;
+
+ d0 = (uint64_t)a[0]*b[0]
+ + 5 * ((uint64_t)a[4]*b[1] + (uint64_t)a[3]*b[2]
+ + (uint64_t)a[2]*b[3] + (uint64_t)a[1]*b[4]);
+ d1 = (uint64_t)a[0]*b[1] + (uint64_t)a[1]*b[0]
+ + 5 * ((uint64_t)a[4]*b[2] + (uint64_t)a[3]*b[3]
+ + (uint64_t)a[2]*b[4]);
+ d2 = (uint64_t)a[0]*b[2] + (uint64_t)a[1]*b[1]
+ + (uint64_t)a[2]*b[0]
+ + 5 * ((uint64_t)a[4]*b[3] + (uint64_t)a[3]*b[4]);
+ d3 = (uint64_t)a[0]*b[3] + (uint64_t)a[1]*b[2]
+ + (uint64_t)a[2]*b[1] + (uint64_t)a[3]*b[0]
+ + 5 * ((uint64_t)a[4]*b[4]);
+ d4 = (uint64_t)a[0]*b[4] + (uint64_t)a[1]*b[3]
+ + (uint64_t)a[2]*b[2] + (uint64_t)a[3]*b[1]
+ + (uint64_t)a[4]*b[0];
+
+ c = d0 >> 26; d0 &= MASK26; d1 += c;
+ c = d1 >> 26; d1 &= MASK26; d2 += c;
+ c = d2 >> 26; d2 &= MASK26; d3 += c;
+ c = d3 >> 26; d3 &= MASK26; d4 += c;
+ c = d4 >> 26; d4 &= MASK26; d0 += c * 5;
+ c = d0 >> 26; d0 &= MASK26; d1 += c;
+
+ out[0] = (uint32_t)d0;
+ out[1] = (uint32_t)d1;
+ out[2] = (uint32_t)d2;
+ out[3] = (uint32_t)d3;
+ out[4] = (uint32_t)d4;
+}
+
+/*
+ * Parse 16 little-endian bytes plus a 'hibit' value (0 or 1) at bit
+ * 128 into 5 26-bit limbs.
+ */
+static inline void blk2limbs(const uint8_t m[16], uint32_t hibit,
+ uint32_t l[5]) {
+ uint32_t t0, t1, t2, t3;
+ memcpy(&t0, m, 4);
+ memcpy(&t1, m + 4, 4);
+ memcpy(&t2, m + 8, 4);
+ memcpy(&t3, m + 12, 4);
+ l[0] = t0 & MASK26;
+ l[1] = ((t0 >> 26) | (t1 << 6)) & MASK26;
+ l[2] = ((t1 >> 20) | (t2 << 12)) & MASK26;
+ l[3] = ((t2 >> 14) | (t3 << 18)) & MASK26;
+ l[4] = (t3 >> 8) | (hibit << 24);
+}
+
+/*
+ * Process one full 16-byte block: h := (h + m) * r mod p, scalar
+ * implementation.
+ */
+static inline void scalar_block(uint32_t h[5], const uint32_t r[5],
+ const uint8_t m[16], uint32_t hibit) {
+ uint32_t blk[5];
+ blk2limbs(m, hibit, blk);
+
+ uint32_t hl[5];
+ hl[0] = h[0] + blk[0];
+ hl[1] = h[1] + blk[1];
+ hl[2] = h[2] + blk[2];
+ hl[3] = h[3] + blk[3];
+ hl[4] = h[4] + blk[4];
+
+ mul_mod_p(hl, r, h);
+}
+
+#if defined(__aarch64__)
+
+/*
+ * 4-block NEON update. Computes:
+ * h := (h + m_0)*r^4 + m_1*r^3 + m_2*r^2 + m_3*r^1 mod p
+ *
+ * where m_0..m_3 are four consecutive 16-byte input blocks (so the
+ * function consumes 64 bytes). The polynomial identity is the
+ * standard Horner expansion of 4 sequential block updates.
+ */
+static void neon4_block(uint32_t h[5],
+ const uint32_t r1[5], const uint32_t r2[5],
+ const uint32_t r3[5], const uint32_t r4[5],
+ const uint8_t m[64]) {
+ /* Limbify the four blocks. */
+ uint32_t b0[5], b1[5], b2[5], b3[5];
+ blk2limbs(m, 1, b0);
+ blk2limbs(m + 16, 1, b1);
+ blk2limbs(m + 32, 1, b2);
+ blk2limbs(m + 48, 1, b3);
+
+ /* Fold the running accumulator into the first block. */
+ b0[0] += h[0]; b0[1] += h[1]; b0[2] += h[2];
+ b0[3] += h[3]; b0[4] += h[4];
+
+ /* Pack messages: mv[i] = (b0[i], b1[i], b2[i], b3[i]). */
+ uint32x4_t mv0 = { b0[0], b1[0], b2[0], b3[0] };
+ uint32x4_t mv1 = { b0[1], b1[1], b2[1], b3[1] };
+ uint32x4_t mv2 = { b0[2], b1[2], b2[2], b3[2] };
+ uint32x4_t mv3 = { b0[3], b1[3], b2[3], b3[3] };
+ uint32x4_t mv4 = { b0[4], b1[4], b2[4], b3[4] };
+
+ /* Pack r-powers: rv[i] = (r^4[i], r^3[i], r^2[i], r^1[i]). */
+ uint32x4_t rv0 = { r4[0], r3[0], r2[0], r1[0] };
+ uint32x4_t rv1 = { r4[1], r3[1], r2[1], r1[1] };
+ uint32x4_t rv2 = { r4[2], r3[2], r2[2], r1[2] };
+ uint32x4_t rv3 = { r4[3], r3[3], r2[3], r1[3] };
+ uint32x4_t rv4 = { r4[4], r3[4], r2[4], r1[4] };
+
+ /* 5 * r-powers, for partial products whose position wraps past
+ * 4 (mod 2^130 - 5 = 5). */
+ uint32x4_t rv1_5 = vaddq_u32(vshlq_n_u32(rv1, 2), rv1);
+ uint32x4_t rv2_5 = vaddq_u32(vshlq_n_u32(rv2, 2), rv2);
+ uint32x4_t rv3_5 = vaddq_u32(vshlq_n_u32(rv3, 2), rv3);
+ uint32x4_t rv4_5 = vaddq_u32(vshlq_n_u32(rv4, 2), rv4);
+
+ /*
+ * Output limb j = sum_i (mv[i] * appropriate r-power for shift j-i,
+ * with 5x multiplier when j - i is negative). Per output we sum
+ * across the 4 lanes ('vaddvq_u64') after collecting both vmull
+ * halves.
+ */
+#define MUL5_LO(m_, r_) vmull_u32(vget_low_u32(m_), vget_low_u32(r_))
+#define MUL5_HI(m_, r_) vmull_high_u32(m_, r_)
+#define MLA_LO(acc, m_, r_) \
+ vmlal_u32(acc, vget_low_u32(m_), vget_low_u32(r_))
+#define MLA_HI(acc, m_, r_) vmlal_high_u32(acc, m_, r_)
+
+ uint64x2_t l0 = MUL5_LO(mv0, rv0);
+ uint64x2_t h0 = MUL5_HI(mv0, rv0);
+ l0 = MLA_LO(l0, mv1, rv4_5); h0 = MLA_HI(h0, mv1, rv4_5);
+ l0 = MLA_LO(l0, mv2, rv3_5); h0 = MLA_HI(h0, mv2, rv3_5);
+ l0 = MLA_LO(l0, mv3, rv2_5); h0 = MLA_HI(h0, mv3, rv2_5);
+ l0 = MLA_LO(l0, mv4, rv1_5); h0 = MLA_HI(h0, mv4, rv1_5);
+ uint64_t d0 = vaddvq_u64(l0) + vaddvq_u64(h0);
+
+ uint64x2_t l1 = MUL5_LO(mv0, rv1);
+ uint64x2_t h1 = MUL5_HI(mv0, rv1);
+ l1 = MLA_LO(l1, mv1, rv0 ); h1 = MLA_HI(h1, mv1, rv0 );
+ l1 = MLA_LO(l1, mv2, rv4_5); h1 = MLA_HI(h1, mv2, rv4_5);
+ l1 = MLA_LO(l1, mv3, rv3_5); h1 = MLA_HI(h1, mv3, rv3_5);
+ l1 = MLA_LO(l1, mv4, rv2_5); h1 = MLA_HI(h1, mv4, rv2_5);
+ uint64_t d1 = vaddvq_u64(l1) + vaddvq_u64(h1);
+
+ uint64x2_t l2 = MUL5_LO(mv0, rv2);
+ uint64x2_t h2 = MUL5_HI(mv0, rv2);
+ l2 = MLA_LO(l2, mv1, rv1 ); h2 = MLA_HI(h2, mv1, rv1 );
+ l2 = MLA_LO(l2, mv2, rv0 ); h2 = MLA_HI(h2, mv2, rv0 );
+ l2 = MLA_LO(l2, mv3, rv4_5); h2 = MLA_HI(h2, mv3, rv4_5);
+ l2 = MLA_LO(l2, mv4, rv3_5); h2 = MLA_HI(h2, mv4, rv3_5);
+ uint64_t d2 = vaddvq_u64(l2) + vaddvq_u64(h2);
+
+ uint64x2_t l3 = MUL5_LO(mv0, rv3);
+ uint64x2_t h3 = MUL5_HI(mv0, rv3);
+ l3 = MLA_LO(l3, mv1, rv2 ); h3 = MLA_HI(h3, mv1, rv2 );
+ l3 = MLA_LO(l3, mv2, rv1 ); h3 = MLA_HI(h3, mv2, rv1 );
+ l3 = MLA_LO(l3, mv3, rv0 ); h3 = MLA_HI(h3, mv3, rv0 );
+ l3 = MLA_LO(l3, mv4, rv4_5); h3 = MLA_HI(h3, mv4, rv4_5);
+ uint64_t d3 = vaddvq_u64(l3) + vaddvq_u64(h3);
+
+ uint64x2_t l4 = MUL5_LO(mv0, rv4);
+ uint64x2_t h4 = MUL5_HI(mv0, rv4);
+ l4 = MLA_LO(l4, mv1, rv3); h4 = MLA_HI(h4, mv1, rv3);
+ l4 = MLA_LO(l4, mv2, rv2); h4 = MLA_HI(h4, mv2, rv2);
+ l4 = MLA_LO(l4, mv3, rv1); h4 = MLA_HI(h4, mv3, rv1);
+ l4 = MLA_LO(l4, mv4, rv0); h4 = MLA_HI(h4, mv4, rv0);
+ uint64_t d4 = vaddvq_u64(l4) + vaddvq_u64(h4);
+
+#undef MUL5_LO
+#undef MUL5_HI
+#undef MLA_LO
+#undef MLA_HI
+
+ /* Carry propagation, same shape as 'mul_mod_p'. */
+ uint64_t c;
+ c = d0 >> 26; d0 &= MASK26; d1 += c;
+ c = d1 >> 26; d1 &= MASK26; d2 += c;
+ c = d2 >> 26; d2 &= MASK26; d3 += c;
+ c = d3 >> 26; d3 &= MASK26; d4 += c;
+ c = d4 >> 26; d4 &= MASK26; d0 += c * 5;
+ c = d0 >> 26; d0 &= MASK26; d1 += c;
+
+ h[0] = (uint32_t)d0;
+ h[1] = (uint32_t)d1;
+ h[2] = (uint32_t)d2;
+ h[3] = (uint32_t)d3;
+ h[4] = (uint32_t)d4;
+}
+
+#elif defined(__x86_64__)
+
+/*
+ * Return 1 if the CPU implements AVX2 and the OS preserves the YMM
+ * registers across context switches, 0 otherwise.
+ */
+static int avx2_supported(void) {
+ unsigned int eax, ebx, ecx, edx, xcr0_lo, xcr0_hi;
+ if (!__get_cpuid(1, &eax, &ebx, &ecx, &edx))
+ return 0;
+ if (!(ecx & bit_OSXSAVE) || !(ecx & bit_AVX))
+ return 0;
+ /* XCR0 bits 1 and 2: XMM and YMM state enabled by the OS */
+ __asm__("xgetbv" : "=a"(xcr0_lo), "=d"(xcr0_hi) : "c"(0));
+ if ((xcr0_lo & 6) != 6)
+ return 0;
+ if (!__get_cpuid_count(7, 0, &eax, &ebx, &ecx, &edx))
+ return 0;
+ return (ebx & bit_AVX2) != 0;
+}
+
+/*
+ * 'avx2_supported', computed once. CPUID is slow, and may trap to
+ * the hypervisor under virtualisation. Racing first calls store the
+ * same value, so relaxed atomics suffice.
+ */
+static int have_avx2(void) {
+ static int cached = -1;
+ int v = __atomic_load_n(&cached, __ATOMIC_RELAXED);
+ if (v < 0) {
+ v = avx2_supported();
+ __atomic_store_n(&cached, v, __ATOMIC_RELAXED);
+ }
+ return v;
+}
+
+/* Sum the four 64-bit lanes of 'v'. */
+__attribute__((target("avx2")))
+static inline uint64_t hsum4(__m256i v) {
+ __m128i s = _mm_add_epi64(_mm256_castsi256_si128(v),
+ _mm256_extracti128_si256(v, 1));
+ s = _mm_add_epi64(s, _mm_unpackhi_epi64(s, s));
+ return (uint64_t)_mm_cvtsi128_si64(s);
+}
+
+/*
+ * AVX2 4-block update, computing the same thing as the NEON kernel:
+ *
+ * h := (h + m_0)*r^4 + m_1*r^3 + m_2*r^2 + m_3*r^1 mod p
+ *
+ * for each successive 64-byte group of 'msg', while 'msglen' allows.
+ * Returns the number of bytes consumed.
+ *
+ * Each '__m256i' holds one limb position across four blocks, one limb
+ * per 64-bit lane. '_mm256_mul_epu32' multiplies the low 32 bits of
+ * each lane into a 64-bit product, so it plays the role of NEON's
+ * 'vmull'/'vmlal' pairs, with a single vector covering all 4 lanes.
+ *
+ * Compiled for AVX2 only; callers must first check 'have_avx2'.
+ */
+__attribute__((target("avx2")))
+static size_t avx2_blocks(uint32_t h[5],
+ const uint32_t r1[5], const uint32_t r2[5],
+ const uint32_t r3[5], const uint32_t r4[5],
+ const uint8_t *msg, size_t msglen) {
+ /* rv[i] = (r^4[i], r^3[i], r^2[i], r^1[i]), lowest lane first. */
+ __m256i rv0 = _mm256_setr_epi64x(r4[0], r3[0], r2[0], r1[0]);
+ __m256i rv1 = _mm256_setr_epi64x(r4[1], r3[1], r2[1], r1[1]);
+ __m256i rv2 = _mm256_setr_epi64x(r4[2], r3[2], r2[2], r1[2]);
+ __m256i rv3 = _mm256_setr_epi64x(r4[3], r3[3], r2[3], r1[3]);
+ __m256i rv4 = _mm256_setr_epi64x(r4[4], r3[4], r2[4], r1[4]);
+
+ /* 5 * r-powers, for partial products whose position wraps past
+ * 4 (mod 2^130 - 5 = 5). */
+ __m256i rv1_5 = _mm256_add_epi64(_mm256_slli_epi64(rv1, 2), rv1);
+ __m256i rv2_5 = _mm256_add_epi64(_mm256_slli_epi64(rv2, 2), rv2);
+ __m256i rv3_5 = _mm256_add_epi64(_mm256_slli_epi64(rv3, 2), rv3);
+ __m256i rv4_5 = _mm256_add_epi64(_mm256_slli_epi64(rv4, 2), rv4);
+
+ size_t pos = 0;
+ while (pos + 64 <= msglen) {
+ const uint8_t *m = msg + pos;
+
+ /* Limbify the four blocks. */
+ uint32_t b0[5], b1[5], b2[5], b3[5];
+ blk2limbs(m, 1, b0);
+ blk2limbs(m + 16, 1, b1);
+ blk2limbs(m + 32, 1, b2);
+ blk2limbs(m + 48, 1, b3);
+
+ /* Fold the running accumulator into the first block. */
+ b0[0] += h[0]; b0[1] += h[1]; b0[2] += h[2];
+ b0[3] += h[3]; b0[4] += h[4];
+
+ /* mv[i] = (b0[i], b1[i], b2[i], b3[i]), lowest lane first. */
+ __m256i mv0 = _mm256_setr_epi64x(b0[0], b1[0], b2[0], b3[0]);
+ __m256i mv1 = _mm256_setr_epi64x(b0[1], b1[1], b2[1], b3[1]);
+ __m256i mv2 = _mm256_setr_epi64x(b0[2], b1[2], b2[2], b3[2]);
+ __m256i mv3 = _mm256_setr_epi64x(b0[3], b1[3], b2[3], b3[3]);
+ __m256i mv4 = _mm256_setr_epi64x(b0[4], b1[4], b2[4], b3[4]);
+
+#define MUL(m_, r_) _mm256_mul_epu32(m_, r_)
+#define MLA(acc, m_, r_) _mm256_add_epi64(acc, _mm256_mul_epu32(m_, r_))
+
+ __m256i t;
+
+ t = MUL(mv0, rv0);
+ t = MLA(t, mv1, rv4_5);
+ t = MLA(t, mv2, rv3_5);
+ t = MLA(t, mv3, rv2_5);
+ t = MLA(t, mv4, rv1_5);
+ uint64_t d0 = hsum4(t);
+
+ t = MUL(mv0, rv1);
+ t = MLA(t, mv1, rv0);
+ t = MLA(t, mv2, rv4_5);
+ t = MLA(t, mv3, rv3_5);
+ t = MLA(t, mv4, rv2_5);
+ uint64_t d1 = hsum4(t);
+
+ t = MUL(mv0, rv2);
+ t = MLA(t, mv1, rv1);
+ t = MLA(t, mv2, rv0);
+ t = MLA(t, mv3, rv4_5);
+ t = MLA(t, mv4, rv3_5);
+ uint64_t d2 = hsum4(t);
+
+ t = MUL(mv0, rv3);
+ t = MLA(t, mv1, rv2);
+ t = MLA(t, mv2, rv1);
+ t = MLA(t, mv3, rv0);
+ t = MLA(t, mv4, rv4_5);
+ uint64_t d3 = hsum4(t);
+
+ t = MUL(mv0, rv4);
+ t = MLA(t, mv1, rv3);
+ t = MLA(t, mv2, rv2);
+ t = MLA(t, mv3, rv1);
+ t = MLA(t, mv4, rv0);
+ uint64_t d4 = hsum4(t);
+
+#undef MUL
+#undef MLA
+
+ /* Carry propagation, same shape as 'mul_mod_p'. */
+ uint64_t c;
+ c = d0 >> 26; d0 &= MASK26; d1 += c;
+ c = d1 >> 26; d1 &= MASK26; d2 += c;
+ c = d2 >> 26; d2 &= MASK26; d3 += c;
+ c = d3 >> 26; d3 &= MASK26; d4 += c;
+ c = d4 >> 26; d4 &= MASK26; d0 += c * 5;
+ c = d0 >> 26; d0 &= MASK26; d1 += c;
+
+ h[0] = (uint32_t)d0;
+ h[1] = (uint32_t)d1;
+ h[2] = (uint32_t)d2;
+ h[3] = (uint32_t)d3;
+ h[4] = (uint32_t)d4;
+
+ pos += 64;
+ }
+ return pos;
+}
+
+#endif
+
+/*
+ * Compute a 16-byte Poly1305 MAC over 'msg' (length 'msglen') using
+ * the 32-byte 'key'. Writes the tag to 'mac_out'.
+ */
+void poly1305_mac_native(const uint8_t key[32], const uint8_t *msg,
+ size_t msglen, uint8_t mac_out[16]) {
+ /* clamp r */
+ uint32_t t0, t1, t2, t3;
+ memcpy(&t0, key, 4);
+ memcpy(&t1, key + 4, 4);
+ memcpy(&t2, key + 8, 4);
+ memcpy(&t3, key + 12, 4);
+ t0 &= 0x0fffffffu;
+ t1 &= 0x0ffffffcu;
+ t2 &= 0x0ffffffcu;
+ t3 &= 0x0ffffffcu;
+
+ uint32_t r[5];
+ r[0] = t0 & MASK26;
+ r[1] = ((t0 >> 26) | (t1 << 6)) & MASK26;
+ r[2] = ((t1 >> 20) | (t2 << 12)) & MASK26;
+ r[3] = ((t2 >> 14) | (t3 << 18)) & MASK26;
+ r[4] = (t3 >> 8);
+
+ uint32_t h[5] = { 0, 0, 0, 0, 0 };
+
+ size_t pos = 0;
+
+#if defined(__aarch64__)
+ /* NEON 4-way path: amortizing the r^2/r^3/r^4 precomputation
+ * (3 scalar mul_mod_p calls) needs about 4 NEON iterations to
+ * break even versus the scalar block loop, so only engage it
+ * when we have at least 16 full blocks (256 bytes). */
+ if (msglen >= 256) {
+ uint32_t r2[5], r3[5], r4[5];
+ mul_mod_p(r, r, r2);
+ mul_mod_p(r2, r, r3);
+ mul_mod_p(r3, r, r4);
+
+ while (pos + 64 <= msglen) {
+ neon4_block(h, r, r2, r3, r4, msg + pos);
+ pos += 64;
+ }
+ }
+#elif defined(__x86_64__)
+ /* AVX2 4-way path, where the CPU supports it. The scalar kernel
+ * is slower here relative to the r^2/r^3/r^4 precomputation than
+ * on aarch64, so two AVX2 iterations (128 bytes) already break
+ * even. */
+ if (msglen >= 128 && have_avx2()) {
+ uint32_t r2[5], r3[5], r4[5];
+ mul_mod_p(r, r, r2);
+ mul_mod_p(r2, r, r3);
+ mul_mod_p(r3, r, r4);
+
+ pos = avx2_blocks(h, r, r2, r3, r4, msg, msglen);
+ }
+#endif
+
+ /* Scalar tail: any remaining full blocks (< 4 of them). */
+ while (pos + 16 <= msglen) {
+ scalar_block(h, r, msg + pos, 1);
+ pos += 16;
+ }
+
+ /* Final partial block (1..15 trailing bytes), if any. */
+ if (pos < msglen) {
+ size_t rem = msglen - pos;
+ uint8_t pad[16] = { 0 };
+ memcpy(pad, msg + pos, rem);
+ pad[rem] = 1;
+ scalar_block(h, r, pad, 0);
+ }
+
+ /* normalize h (mul_mod_p may leave a small excess in h[1]) */
+ {
+ uint32_t c;
+ c = h[1] >> 26; h[1] &= MASK26; h[2] += c;
+ c = h[2] >> 26; h[2] &= MASK26; h[3] += c;
+ c = h[3] >> 26; h[3] &= MASK26; h[4] += c;
+ c = h[4] >> 26; h[4] &= MASK26; h[0] += c * 5;
+ c = h[0] >> 26; h[0] &= MASK26; h[1] += c;
+ }
+
+ /* full reduction to [0, p) via constant-time conditional
+ * subtraction of p */
+ uint32_t g[5];
+ uint32_t c = 5;
+ g[0] = h[0] + c; c = g[0] >> 26; g[0] &= MASK26;
+ g[1] = h[1] + c; c = g[1] >> 26; g[1] &= MASK26;
+ g[2] = h[2] + c; c = g[2] >> 26; g[2] &= MASK26;
+ g[3] = h[3] + c; c = g[3] >> 26; g[3] &= MASK26;
+ g[4] = h[4] + c;
+ uint32_t carry = g[4] >> 26;
+ g[4] &= MASK26;
+ uint32_t mask = (uint32_t)0 - carry;
+ h[0] = (h[0] & ~mask) | (g[0] & mask);
+ h[1] = (h[1] & ~mask) | (g[1] & mask);
+ h[2] = (h[2] & ~mask) | (g[2] & mask);
+ h[3] = (h[3] & ~mask) | (g[3] & mask);
+ h[4] = (h[4] & ~mask) | (g[4] & mask);
+
+ /* repack 5x 26-bit limbs into 4x 32-bit limbs */
+ uint32_t h0 = h[0] | (h[1] << 26);
+ uint32_t h1 = (h[1] >> 6) | (h[2] << 20);
+ uint32_t h2 = (h[2] >> 12) | (h[3] << 14);
+ uint32_t h3 = (h[3] >> 18) | (h[4] << 8);
+
+ /* add s (high 16 bytes of key), mod 2^128 */
+ uint32_t s0, s1, s2, s3;
+ memcpy(&s0, key + 16, 4);
+ memcpy(&s1, key + 20, 4);
+ memcpy(&s2, key + 24, 4);
+ memcpy(&s3, key + 28, 4);
+
+ uint64_t a0 = (uint64_t)h0 + s0;
+ uint64_t a1 = (uint64_t)h1 + s1 + (a0 >> 32);
+ uint64_t a2 = (uint64_t)h2 + s2 + (a1 >> 32);
+ uint64_t a3 = (uint64_t)h3 + s3 + (a2 >> 32);
+
+ uint32_t o0 = (uint32_t)a0;
+ uint32_t o1 = (uint32_t)a1;
+ uint32_t o2 = (uint32_t)a2;
+ uint32_t o3 = (uint32_t)a3;
+
+ memcpy(mac_out + 0, &o0, 4);
+ memcpy(mac_out + 4, &o1, 4);
+ memcpy(mac_out + 8, &o2, 4);
+ memcpy(mac_out + 12, &o3, 4);
+}
+
+int poly1305_native_available(void) {
+ return 1;
+}
+
+#else
+
+void poly1305_mac_native(const uint8_t *key, const uint8_t *msg,
+ size_t msglen, uint8_t *mac_out) {
+ (void)key; (void)msg; (void)msglen; (void)mac_out;
+}
+
+int poly1305_native_available(void) {
+ return 0;
+}
+
+#endif
diff --git a/lib/Crypto/MAC/Poly1305.hs b/lib/Crypto/MAC/Poly1305.hs
@@ -26,7 +26,7 @@ module Crypto.MAC.Poly1305 (
, _roll16
) where
-import qualified Crypto.MAC.Poly1305.Arm as Arm
+import qualified Crypto.MAC.Poly1305.Native as Native
import Data.Barrier (barrier)
import qualified Data.Bits as B
import qualified Data.ByteString as BS
@@ -189,8 +189,8 @@ mac
-> Maybe MAC -- ^ 128-bit message authentication code
mac key@(BI.PS _ _ kl) msg
| kl /= 32 = Nothing
- | Arm.poly1305_arm_available =
- pure $! MAC (Arm.mac key msg)
+ | Native.poly1305_native_available =
+ pure $! MAC (Native.mac key msg)
| otherwise =
let (clamp . _roll16 -> r, _roll16 -> s) = BS.splitAt 16 key
in pure $! (MAC (_poly1305_loop r s msg))
diff --git a/lib/Crypto/MAC/Poly1305/Arm.hs b/lib/Crypto/MAC/Poly1305/Arm.hs
@@ -1,58 +0,0 @@
-{-# OPTIONS_HADDOCK hide #-}
-{-# LANGUAGE BangPatterns #-}
-
--- |
--- Module: Crypto.MAC.Poly1305.Arm
--- Copyright: (c) 2025 Jared Tobin
--- License: MIT
--- Maintainer: Jared Tobin <jared@ppad.tech>
---
--- ARM acceleration for the Poly1305 MAC.
-
-module Crypto.MAC.Poly1305.Arm (
- poly1305_arm_available
- , mac
- ) where
-
-import qualified Data.ByteString as BS
-import qualified Data.ByteString.Internal as BI
-import Data.Word (Word8)
-import Foreign.C.Types (CInt(..), CSize(..))
-import Foreign.ForeignPtr (withForeignPtr)
-import Foreign.Ptr (Ptr, plusPtr)
-import System.IO.Unsafe (unsafeDupablePerformIO)
-
--- ffi ------------------------------------------------------------------------
-
-foreign import ccall unsafe "poly1305_mac_arm"
- c_poly1305_mac
- :: Ptr Word8 -> Ptr Word8 -> CSize -> Ptr Word8 -> IO ()
-
-foreign import ccall unsafe "poly1305_arm_available"
- c_poly1305_arm_available :: IO CInt
-
--- utilities ------------------------------------------------------------------
-
-fi :: (Integral a, Num b) => a -> b
-fi = fromIntegral
-{-# INLINE fi #-}
-
--- api ------------------------------------------------------------------------
-
--- | Are ARM extensions available?
-poly1305_arm_available :: Bool
-poly1305_arm_available =
- unsafeDupablePerformIO c_poly1305_arm_available /= 0
-{-# NOINLINE poly1305_arm_available #-}
-
--- | Compute a Poly1305 MAC over the message using the given (already-
--- validated 32-byte) key.
-mac :: BS.ByteString -> BS.ByteString -> BS.ByteString
-mac (BI.PS kfp koff _) (BI.PS mfp moff mlen) =
- BI.unsafeCreate 16 $ \dst ->
- withForeignPtr kfp $ \kp0 ->
- withForeignPtr mfp $ \mp0 ->
- c_poly1305_mac (kp0 `plusPtr` koff)
- (mp0 `plusPtr` moff)
- (fi mlen)
- dst
diff --git a/lib/Crypto/MAC/Poly1305/Native.hs b/lib/Crypto/MAC/Poly1305/Native.hs
@@ -0,0 +1,59 @@
+{-# OPTIONS_HADDOCK hide #-}
+{-# LANGUAGE BangPatterns #-}
+
+-- |
+-- Module: Crypto.MAC.Poly1305.Native
+-- Copyright: (c) 2025 Jared Tobin
+-- License: MIT
+-- Maintainer: Jared Tobin <jared@ppad.tech>
+--
+-- A native C kernel for the Poly1305 MAC on aarch64 and x86_64, with
+-- ARM NEON or x86 AVX2 acceleration.
+
+module Crypto.MAC.Poly1305.Native (
+ poly1305_native_available
+ , mac
+ ) where
+
+import qualified Data.ByteString as BS
+import qualified Data.ByteString.Internal as BI
+import Data.Word (Word8)
+import Foreign.C.Types (CInt(..), CSize(..))
+import Foreign.ForeignPtr (withForeignPtr)
+import Foreign.Ptr (Ptr, plusPtr)
+import System.IO.Unsafe (unsafeDupablePerformIO)
+
+-- ffi ------------------------------------------------------------------------
+
+foreign import ccall unsafe "poly1305_mac_native"
+ c_poly1305_mac
+ :: Ptr Word8 -> Ptr Word8 -> CSize -> Ptr Word8 -> IO ()
+
+foreign import ccall unsafe "poly1305_native_available"
+ c_poly1305_native_available :: IO CInt
+
+-- utilities ------------------------------------------------------------------
+
+fi :: (Integral a, Num b) => a -> b
+fi = fromIntegral
+{-# INLINE fi #-}
+
+-- api ------------------------------------------------------------------------
+
+-- | Is the native kernel available?
+poly1305_native_available :: Bool
+poly1305_native_available =
+ unsafeDupablePerformIO c_poly1305_native_available /= 0
+{-# NOINLINE poly1305_native_available #-}
+
+-- | Compute a Poly1305 MAC over the message using the given (already-
+-- validated 32-byte) key.
+mac :: BS.ByteString -> BS.ByteString -> BS.ByteString
+mac (BI.PS kfp koff _) (BI.PS mfp moff mlen) =
+ BI.unsafeCreate 16 $ \dst ->
+ withForeignPtr kfp $ \kp0 ->
+ withForeignPtr mfp $ \mp0 ->
+ c_poly1305_mac (kp0 `plusPtr` koff)
+ (mp0 `plusPtr` moff)
+ (fi mlen)
+ dst
diff --git a/ppad-poly1305.cabal b/ppad-poly1305.cabal
@@ -44,14 +44,14 @@ library
exposed-modules:
Crypto.MAC.Poly1305
other-modules:
- Crypto.MAC.Poly1305.Arm
+ Crypto.MAC.Poly1305.Native
Data.Barrier
build-depends:
base >= 4.9 && < 5
, bytestring >= 0.9 && < 0.13
, ppad-fixed >= 0.2.1 && < 0.3
c-sources:
- cbits/poly1305_arm.c
+ cbits/poly1305_native.c
if arch(aarch64)
cc-options: -march=armv8-a
if flag(sanitize)
diff --git a/test/Main.hs b/test/Main.hs
@@ -6,7 +6,10 @@
module Main where
import qualified Crypto.MAC.Poly1305 as Poly1305
+import Data.Bits ((.&.))
+import qualified Data.ByteString as BS
import qualified Data.ByteString.Base16 as B16
+import Data.Foldable (for_)
import Test.Tasty
import qualified Test.Tasty.HUnit as H
@@ -24,6 +27,7 @@ main = defaultMain $ testGroup "ppad-poly1305" [
, mac9
, mac10
, mac11
+ , dispatch
]
mac :: TestTree
@@ -172,3 +176,29 @@ mac11 = H.testCase "mac (A.3 #11)" $ do
out = Poly1305._poly1305_loop r s msg
H.assertEqual mempty tag out
+-- dispatching vs pure
+
+-- the dispatching 'mac' agrees with the pure Haskell loop for every
+-- message length from 0 to 640 bytes, covering partial blocks, short
+-- tails, and any vectorised paths
+dispatch :: TestTree
+dispatch = H.testCase "mac (dispatch vs pure)" $
+ for_ keys $ \key ->
+ for_ msgs $ \m ->
+ for_ [0..640] $ \l -> do
+ let msg = BS.take l m
+ Just out = Poly1305.mac key msg
+ H.assertEqual ("length " <> show l) (Poly1305.MAC (ref key msg)) out
+ where
+ Just clamp = B16.decode "ffffff0ffcffff0ffcffff0ffcffff0f"
+ Just key0 = B16.decode
+ "85d6be7857556d337f4452fe42d506a80103808afb0db2fd4abff6af4149f51b"
+ keys = [key0, BS.replicate 32 0xff]
+ msgs = [ BS.pack [fromIntegral (i * 7 + 3) | i <- [0..639 :: Int]]
+ , BS.replicate 640 0xff
+ ]
+ ref key msg =
+ let (k0, k1) = BS.splitAt 16 key
+ r = Poly1305._roll16 (BS.packZipWith (.&.) clamp k0)
+ s = Poly1305._roll16 k1
+ in Poly1305._poly1305_loop r s msg
