aere-research/pq-stark/poseidon2_babybear_reference.py
Aere Network 6cb0140fae Republished from a clean root: the compiled artifact is gone from history, and the local line of work joins the sanitized public line
The public history carried kat/__pycache__/mlkem768_reference.cpython-314.pyc,
a compiled Python artifact embedding the operator's absolute local path. Text
secret scanners do not read compiled binaries, which is exactly how it slipped
through, and removing it from the tip would have left it reachable through the
old root commits. So this repository is republished from a single clean root.

This root also carries, from the previously unpublished line of work:
- corrected LICENSE year, LICENSING.md, VERIFY-POLICY.md, and
  CITATIONS-UNRESOLVED.md remeasured 2026-08-11 (101 paths, README aligned)
- O-018: run_consensus_verification.py ran 19 of 29 models and reported PASS;
  it now runs all 29, and computemarket_smt.py gains resolveByTimeout /
  reclaimUnsettled cases plus a negative control
- O-006: the word 'audited' removed from next to Bouncy Castle, twice, after a
  concurrent edit resurrected it
- O-014: prior art named and dated - Algorand's native falcon_verify shipped
  about ten months before AERE's precompiles; the primacy claim is withdrawn
  where it was implied
- bench/ scripts parametrized so they actually run for an outsider (the
  earlier textual sanitization left $STAGING unexpanded inside Python strings)
- AIP-2/AIP-3 errata with measured figures, spec remeasurements at 2026-08-01,
  and the spec-zk-stack retractions (owner is an operational key, not the
  Foundation; 'maximally sound' withdrawn; aggregator V1 deprecated)
The redacted bench-host environment files from the sanitized line are kept
exactly as published; the unredacted local variants are not carried.
2026-08-15 13:52:14 +03:00

407 lines
18 KiB
Python

#!/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)