#!/usr/bin/env python3 # Independent reference for the Poseidon2 permutation over BabyBear at width 16, component (b) of the # PQ STARK-verify precompile 0x0AE8 (direct SP1/Plonky3 inner FRI/STARK verify). Poseidon2 is the # hash Plonky3/SP1 uses BOTH as the Merkle/MMCS compression AND as the Fiat-Shamir duplex-sponge # challenger. See docs/AERE-STARK-VERIFIER-PORT-SPEC.md section 3. # # ============================ HONEST SCOPE (read first) ============================ # This module implements the Poseidon2 permutation and its constants are now CONFIRMED-FROM-SOURCE # against the pinned Plonky3 revision, and a real known-answer test PASSES (see below). What that # does and does NOT mean: # - CONFIRMED (2026-07-19): the 141 round constants (128 external + 13 internal), the internal # diffusion diagonal, the external MDS-light M4, the round counts (ROUNDS_F=8, ROUNDS_P=13), the # x^7 S-box, and the layer ORDER all match p3-baby-bear / p3-poseidon2 0.4.3-succinct (the exact # crates pinned in aerenew Cargo.lock, checksums d69e6e9a.../52298637...). Confirmation method: # the pinned crates were executed via cargo (a byte-identical build, lockfile checksums matched) # to emit the constants and three known-answer permutation vectors; this reference reproduces all # three EXACTLY. See conformance_kats() / test_poseidon2_babybear.py and spec-poseidon2-constants. # - MONTGOMERY SUBTLETY (the trap the port spec warns about, now RESOLVED): p3-baby-bear's # DiffusionMatrixBabyBear internal layer, as a map on CANONICAL vectors, is R^{-1} * (J + diag(D)) # where R = 2^32 mod p and R^{-1} = 943718400. That is, the whole internal layer output is scaled # by the inverse Montgomery factor (verified by recovering the exact 16x16 canonical matrix from # the pinned library via basis-vector probes: off-diagonal = 943718400 = R^{-1}, and diagonal # I[i][i] = R^{-1}*(1 + D[i]) to the last bit). A prior pass used the textbook internal layer # state[i]*D[i] + sum (missing the R^{-1}), which is self-consistent but DISAGREES with the real # prover; that is exactly the "two wrong copies agree" failure mode. It is fixed here. # - The x^7 S-box degree is number-theoretically CONFIRMED (smallest d>1 with gcd(d, p-1)=1 for # BabyBear; proven in the harness), so x^7 is a bijection on F_p. # # What this module still does NOT do: it does not verify any STARK proof, and it is deliberately NOT # wired to make the top-level 0x0AE8 accept anything. The permutation being confirmed conformant is # necessary but not sufficient: the duplex-sponge challenger, the FRI folding/consistency arithmetic, # the MMCS/Merkle openings, and the SP1 recursion-AIR constraint evaluation are still un-ported, so # the precompile stays fail-closed (returns EMPTY for every input). See the README and port spec. # # Source (pinned): p3-poseidon2 0.4.3-succinct (checksum # 522986377b2164c5f94f2dae88e0e0a3d169cc6239202ef4aeb4322d60feffd0) and p3-baby-bear 0.4.3-succinct # (checksum d69e6e9af4eaaaa60f7bb9f0e0f73ebcbaefe7e00974d97ad0fa542d6a4f0890), crates.io registry, # the Plonky3 revision pinned in aerenew/zk-circuits/*/Cargo.lock and rollup-evm-validity/*/Cargo.lock. # Constants generated by Xoroshiro128Plus::seed_from_u64(1) via Poseidon2::new_from_rng_128 with the # Poseidon2ExternalMatrixGeneral external layer and DiffusionMatrixBabyBear internal layer. P = 2013265921 # BabyBear prime 2^31 - 2^27 + 1 = 15*2^27 + 1 = 0x78000001 # Montgomery inverse factor. p3-baby-bear stores field elements in Montgomery form with R = 2^32; its # internal diffusion layer, as a CANONICAL linear map, carries a factor of R^{-1} = (2^32)^{-1} mod p. R = (1 << 32) % P # 268435454 R_INV = pow(R, -1, P) # 943718400 (verified against the pinned library's extracted matrix) WIDTH = 16 # Poseidon2 state width t (Plonky3 default for the SP1 recursion config). CONFIRMED. SBOX_DEGREE = 7 # x^7: smallest d>1 with gcd(d, p-1)=1 (p-1 = 2^27*3*5). CONFIRMED. ROUNDS_F = 8 # external (full) rounds, split 4 initial + 4 terminal. CONFIRMED (round_numbers). ROUNDS_P = 13 # internal (partial) rounds for BabyBear width 16. CONFIRMED (round_numbers). # External 4x4 MDS matrix M4 (Poseidon2 paper). The width-16 external layer M_E is the MDS-light # block construction built from M4 (diagonal blocks 2*M4, off-diagonal blocks M4). CONFIRMED: the # full 16x16 external matrix recovered from Plonky3's Poseidon2ExternalMatrixGeneral equals this. M4 = [ [2, 3, 1, 1], [1, 2, 3, 1], [1, 1, 2, 3], [3, 1, 1, 2], ] # Internal diffusion diagonal diag(M_I - I) = D: the internal linear layer is # M_I = R^{-1} * (J + diag(D)) where J is the all-ones matrix. D[0] = -2 (stored as p-2); the rest are # small powers of two. CONFIRMED against p3-baby-bear 0.4.3-succinct # POSEIDON2_INTERNAL_MATRIX_DIAG_16_BABYBEAR_MONTY (canonical form). INTERNAL_DIAG_M1_16 = [ P - 2, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 32768, ] # ============================ CONFIRMED round constants ============================ # The 141 Poseidon2-BabyBear width-16 round constants (canonical), generated by # Xoroshiro128Plus::seed_from_u64(1) via new_from_rng_128 and extracted from the pinned crates. # External: ROUNDS_F rounds x WIDTH lanes (added to all lanes). Internal: ROUNDS_P (added to lane 0). RC_EXTERNAL = [ [ 1321363468, 285374923, 858595076, 131742120, 550898981, 109281027, 1548327248, 299186948, 1198120888, 1302311359, 568137078, 1484856917, 1301979945, 725688886, 941758026, 323341913, ], [ 1049323172, 822409348, 1406080127, 1279024384, 214862539, 904628921, 1320747287, 11578228, 1036373712, 1474430466, 1430509860, 111174484, 1124450171, 85382027, 679880882, 243277213, ], [ 1338495990, 1523013347, 1841068573, 578194469, 47683837, 1790441672, 1628061601, 1716216090, 1635810049, 1115145248, 1117524270, 678640014, 1962751651, 1367401392, 11688709, 1950824358, ], [ 528649031, 1937116923, 1460949223, 1193074357, 1221801411, 1183923117, 433505619, 1928933309, 505759755, 285671663, 1047265910, 909281502, 1258966486, 864761693, 307024510, 504858517, ], [ 1467478033, 1754565867, 432187324, 1452390672, 881974300, 550050336, 1447309270, 939419487, 1783112406, 1166910332, 107514714, 580516863, 2003318760, 854475946, 934896823, 994783668, ], [ 1841107561, 438269126, 1550523825, 913322122, 600932628, 583000098, 1262690949, 105797869, 277542016, 170491952, 365854467, 1479645308, 1457660602, 1635879552, 499155053, 741227047, ], [ 651389942, 464828001, 89696107, 360044673, 230330371, 1773129416, 1380150763, 745014723, 793475694, 1361274828, 1443741698, 51616650, 731414218, 1087554954, 1273943885, 311581717, ], [ 702702762, 1473247301, 132108357, 1348260424, 476775430, 1438949459, 2434448, 1349232398, 1954471898, 1762138591, 1271221795, 1593266476, 864488771, 139147729, 1053373910, 422842363, ], ] RC_INTERNAL = [ 402771160, 320708227, 1122772462, 100431997, 202594011, 1226485372, 1088619034, 64118538, 109828860, 724723599, 1662837151, 797753907, 1075635743, ] # ============================ base-field helpers ============================ def mul(a, b): return (a * b) % P def add(a, b): return (a + b) % P def sbox_mono(x): """x^7 in F_p (the BabyBear Poseidon2 S-box).""" x = x % P x2 = mul(x, x) x4 = mul(x2, x2) return mul(x4, mul(x2, x)) # x^4 * x^2 * x = x^7 def _sbox_inverse_exponent(): """d^{-1} mod (p-1) where d = 7; used by the inverse permutation for the round-trip bijection test. Exists because gcd(7, p-1) = 1.""" return pow(SBOX_DEGREE, -1, P - 1) _SBOX_INV_E = _sbox_inverse_exponent() def sbox_mono_inv(y): return pow(y % P, _SBOX_INV_E, P) # ============================ linear layers ============================ def _m4_apply(t): """Apply the 4x4 matrix M4 to a length-4 vector.""" return [ (M4[row][0] * t[0] + M4[row][1] * t[1] + M4[row][2] * t[2] + M4[row][3] * t[3]) % P for row in range(4) ] def external_layer(state): """Poseidon2 external (MDS-light) linear layer M_E for width 16 = 4 blocks of 4. Two steps (Plonky3's construction): (1) apply M4 to each block of 4; (2) add the across-block column sums. Equivalently M_E is the block matrix with diagonal blocks 2*M4 and off-diagonal blocks M4.""" blocks = [_m4_apply(state[4 * b:4 * b + 4]) for b in range(4)] col_sum = [(blocks[0][j] + blocks[1][j] + blocks[2][j] + blocks[3][j]) % P for j in range(4)] out = [0] * WIDTH for b in range(4): for j in range(4): out[4 * b + j] = (blocks[b][j] + col_sum[j]) % P return out def internal_layer(state): """Poseidon2 internal diffusion layer, matching p3-baby-bear DiffusionMatrixBabyBear as a canonical map: M_I = R^{-1} * (J + diag(D)), i.e. out[i] = R_INV * (sum + D[i]*state[i]). The R_INV = 943718400 factor is the Montgomery-form artifact of the pinned library (verified by recovering the exact 16x16 canonical matrix); omitting it silently disagrees with the prover.""" s = sum(state) % P return [(R_INV * (s + INTERNAL_DIAG_M1_16[i] * state[i])) % P for i in range(WIDTH)] # ---- explicit matrices (for invertibility / inverse-permutation tests only) ---- def external_matrix(): """The dense 16x16 form of M_E (diagonal blocks 2*M4, off-diagonal blocks M4). Used by the harness to check invertibility and to invert the layer for the bijection round-trip test.""" m = [[0] * WIDTH for _ in range(WIDTH)] for bi in range(4): for bj in range(4): coef = 2 if bi == bj else 1 for r in range(4): for c in range(4): m[4 * bi + r][4 * bj + c] = (coef * M4[r][c]) % P return m def internal_matrix(): """The dense 16x16 form of M_I = R^{-1} * (J + diag(D)); off-diagonal = R_INV, diagonal = R_INV*(1 + D[i]). Matches the matrix recovered from p3-baby-bear DiffusionMatrixBabyBear.""" m = [[R_INV] * WIDTH for _ in range(WIDTH)] for i in range(WIDTH): m[i][i] = (R_INV * (1 + INTERNAL_DIAG_M1_16[i])) % P return m def mat_vec(m, v): return [sum(m[i][j] * v[j] for j in range(WIDTH)) % P for i in range(WIDTH)] def mat_inverse(m): """Inverse of a WIDTH x WIDTH matrix over F_p via Gauss-Jordan. Raises if singular.""" n = WIDTH a = [[m[i][j] % P for j in range(n)] + [1 if j == i else 0 for j in range(n)] for i in range(n)] for col in range(n): piv = next((r for r in range(col, n) if a[r][col] % P != 0), None) if piv is None: raise ValueError("singular matrix (not invertible over F_p)") a[col], a[piv] = a[piv], a[col] inv_p = pow(a[col][col], -1, P) a[col] = [(x * inv_p) % P for x in a[col]] for r in range(n): if r != col and a[r][col] % P != 0: f = a[r][col] a[r] = [(a[r][k] - f * a[col][k]) % P for k in range(2 * n)] return [[a[i][j + n] for j in range(n)] for i in range(n)] # ============================ the permutation ============================ def permute(state): """Poseidon2 permutation over BabyBear, width 16. Structure (Plonky3 order): 0. initial external linear layer (M_E once), 1. ROUNDS_F/2 initial external rounds: add RC to all lanes, S-box all lanes, M_E, 2. ROUNDS_P internal rounds: add RC to lane 0, S-box lane 0, M_I, 3. ROUNDS_F/2 terminal external rounds: add RC to all lanes, S-box all lanes, M_E. Deterministic; a bijection (each layer invertible / S-box a permutation of F_p). CONFIRMED conformant to p3-baby-bear 0.4.3-succinct via known-answer vectors (conformance_kats).""" if len(state) != WIDTH: raise ValueError("state must have WIDTH elements") s = [x % P for x in state] s = external_layer(s) # initial linear layer half = ROUNDS_F // 2 for r in range(half): # initial external rounds s = [add(s[i], RC_EXTERNAL[r][i]) for i in range(WIDTH)] s = [sbox_mono(x) for x in s] s = external_layer(s) for r in range(ROUNDS_P): # internal rounds s[0] = add(s[0], RC_INTERNAL[r]) s[0] = sbox_mono(s[0]) s = internal_layer(s) for r in range(half, ROUNDS_F): # terminal external rounds s = [add(s[i], RC_EXTERNAL[r][i]) for i in range(WIDTH)] s = [sbox_mono(x) for x in s] s = external_layer(s) return s def permute_inverse(state): """Inverse of `permute`, used ONLY to prove the permutation is a genuine bijection (round-trip). Not part of any verifier path.""" if len(state) != WIDTH: raise ValueError("state must have WIDTH elements") s = [x % P for x in state] ext_inv = mat_inverse(external_matrix()) int_inv = mat_inverse(internal_matrix()) half = ROUNDS_F // 2 for r in range(ROUNDS_F - 1, half - 1, -1): # undo terminal external rounds s = mat_vec(ext_inv, s) s = [sbox_mono_inv(x) for x in s] s = [(s[i] - RC_EXTERNAL[r][i]) % P for i in range(WIDTH)] for r in range(ROUNDS_P - 1, -1, -1): # undo internal rounds s = mat_vec(int_inv, s) s[0] = sbox_mono_inv(s[0]) s[0] = (s[0] - RC_INTERNAL[r]) % P for r in range(half - 1, -1, -1): # undo initial external rounds s = mat_vec(ext_inv, s) s = [sbox_mono_inv(x) for x in s] s = [(s[i] - RC_EXTERNAL[r][i]) % P for i in range(WIDTH)] s = mat_vec(ext_inv, s) # undo the initial linear layer return s # ============================ 2-to-1 compression (structural placeholder) ============================ def compress_2to1(left8, right8): """A width-16 -> 8 truncated-permutation compression: absorb l||r (16 elems), permute, take the first 8 output elements. This is the COMMON Plonky3 TruncatedPermutation shape; the EXACT SP1 Merkle compression convention (padding/rate, which 8 lanes) is [VERIFY] for component (c).""" if len(left8) != 8 or len(right8) != 8: raise ValueError("compress inputs must be length 8") return permute(list(left8) + list(right8))[:8] # ============================ CONFIRMED conformance KATs (from the pinned library) ============================ # Real known-answer vectors emitted by executing p3-baby-bear / p3-poseidon2 0.4.3-succinct (the exact # pinned crates; lockfile checksums matched). This reference reproduces all three EXACTLY, which is # what makes the constants + structure CONFIRMED-conformant (not merely self-consistent). CONFORMANCE_KATS = [ { "name": "zeros", "in": [0] * 16, "out": [1787823396, 953829438, 89382455, 347481625, 1754527224, 916217775, 1056029082, 410644796, 1169123478, 1854704276, 1195829987, 1485264906, 1824644035, 1948268315, 847945433, 190591038], }, { "name": "iota", "in": list(range(16)), "out": [157639285, 1851003038, 1852457045, 1920360618, 779990819, 1080011039, 585017685, 1093051731, 249426030, 967262243, 623744062, 280332881, 1995600430, 1751988435, 317724737, 1895035071], }, { "name": "testvec", "in": [894848333, 1437655012, 1200606629, 1690012884, 71131202, 1749206695, 1717947831, 120589055, 19776022, 42382981, 1831865506, 724844064, 171220207, 1299207443, 227047920, 1783754913], "out": [512585766, 975869435, 1921378527, 1238606951, 899635794, 132650430, 1426417547, 1734425242, 57415409, 67173027, 1535042492, 1318033394, 1070659233, 17258943, 856719028, 1500534995], }, ] def conformance_kats(): """Return (passed, total) after running the pinned-library known-answer vectors.""" passed = 0 for kat in CONFORMANCE_KATS: if permute(list(kat["in"])) == kat["out"]: passed += 1 return passed, len(CONFORMANCE_KATS) # ============================ shared cross-language vector set ============================ # A fixed set of input states all three languages (Python, Node, Java) permute; the harness asserts # byte-identical outputs and that each inverse round-trips. Combined with CONFORMANCE_KATS above, this # validates both the arithmetic (cross-language) AND conformance to Plonky3 (known-answer vectors). def _input_states(): states = [ [0] * WIDTH, list(range(WIDTH)), # [0,1,2,...,15] [P - 1] * WIDTH, [(i * 2654435761) % P for i in range(WIDTH)], [11, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1], ] # a few pseudo-random states from a simple deterministic LCG (identical in every language) x = 123456789 for _ in range(8): st = [] for _ in range(WIDTH): x = (1103515245 * x + 12345) % 2147483648 st.append(x % P) states.append(st) return states def shared_vectors(): rc_ext_flat = [RC_EXTERNAL[r][i] for r in range(ROUNDS_F) for i in range(WIDTH)] states = _input_states() return { "constants": { "P": P, "rInv": R_INV, "width": WIDTH, "sboxDegree": SBOX_DEGREE, "roundsF": ROUNDS_F, "roundsP": ROUNDS_P, "internalDiagM1": INTERNAL_DIAG_M1_16, "m4": [M4[r][c] for r in range(4) for c in range(4)], "rcExternalFlat": rc_ext_flat, "rcInternal": RC_INTERNAL, }, "permute": [ {"in": st, "out": permute(st), "roundtrip": permute_inverse(permute(st)) == [x % P for x in st]} for st in states ], "conformanceKats": [ {"name": k["name"], "in": k["in"], "out": permute(list(k["in"])), "match": permute(list(k["in"])) == k["out"]} for k in CONFORMANCE_KATS ], } if __name__ == "__main__": import json import sys json.dump(shared_vectors(), sys.stdout)