import pandas as pd
from IPython.display import Markdown
"""ML-DSA parameter sets, as given in Table 1 of FIPS 204."""
from dataclasses import dataclass
@dataclass
class MLDSAParams:
name: str
n: int # ring degree, always 256
q: int # modulus
tau: int # number of +-1's in the challenge poly c
gamma1: int # y coefficient range
gamma2: int # low-bits rounding range
k: int # rows of A (dimension of w, s2, t)
l: int # columns of A (dimension of s1, z)
eta: int # secret key coefficient range
d: int # dropped bits of t
omega: int # maximum number of hints
@property
def beta(self) -> int:
return self.tau * self.eta
# FIPS-204 Table 1
PARAMS = [
MLDSAParams("ML-DSA-44", n=256, q=8380417, tau=39, gamma1=2**17,
gamma2=(8380417 - 1) // 88, k=4, l=4, eta=2, d=13, omega=80),
MLDSAParams("ML-DSA-65", n=256, q=8380417, tau=49, gamma1=2**19,
gamma2=(8380417 - 1) // 32, k=6, l=5, eta=4, d=13, omega=55),
MLDSAParams("ML-DSA-87", n=256, q=8380417, tau=60, gamma1=2**19,
gamma2=(8380417 - 1) // 32, k=8, l=7, eta=2, d=13, omega=75),
]
"""P_hint for ML-DSA: probability that a signing attempt is rejected on line 28
of FIPS 204 Algorithm 7 because it produced more than omega hints.
Monte Carlo. Per-coefficient hint probability is conditioned on line 23
having passed (|r0| < gamma2 - beta). a = |(c*t0)_i| alone does not determine
whether the hint fires: the count is a-beta-1 if (c*t0)_i = +a, a-beta if
(c*t0)_i = -a (the window (-gamma2, gamma2] is open at one end, closed at the
other, so the two signs are not interchangeable). Since t0 is symmetric about
zero, the sign is a fair coin flip given a, so the per-coefficient probability
is the average of the two:
P(hint | a) = max(0, a - beta - 1/2) / (2*(gamma2 - beta) - 1), a = |(c*t0)_i|
"""
import math
import numpy as np
# No import of params.py here: this module has no import-time dependency
# on it, so it stays executable whether imported as a package module or
# spliced inline into a notebook cell. run_simulation() below takes the
# parameter list as an explicit argument instead of defaulting to it.
PSI = np.exp(1j * np.pi * np.arange(256) / 256)
IPSI = np.conj(PSI)
def p_hint(k, tau, d, gamma2, beta, omega, trials, batch=2000, seed=19820116):
"""Monte Carlo estimate of P(N > omega), N = number of hints that fire."""
rng = np.random.default_rng(seed)
half = 2 ** (d - 1)
denom = 2 * (gamma2 - beta) - 1
total_hints = rejected = done = ct0_rej = 0
while done < trials:
b = min(batch, trials - done)
# Step 1: sample c
c = np.zeros((b, 256))
# Spread tau's among 256 positions
pos = np.argsort(rng.random((b, 256)), axis=1)[:, :tau]
# Map {0,1} to {-1,+1} and put in chosen positions.
np.put_along_axis(c, pos,
(rng.integers(0, 2, size=(b, tau)) * 2 - 1).astype(float),
axis=1)
# Step 2: sample t0 uniform on (-2^(d-1), 2^(d-1)];
t0 = rng.integers(-half + 1, half + 1, size=(b, k, 256)).astype(float)
# Step 3: compute c*t0 in Z[X]/(X^256+1) via negacyclic FFT
ct0 = np.real(np.fft.ifft(np.fft.fft(t0 * PSI, axis=2)
* np.fft.fft(c * PSI, axis=1)[:, None, :],
axis=2) * IPSI)
# Number of rejected ct0. We expect 0
ct0_rej += int((np.abs(np.rint(ct0)).max(axis=(1, 2)) >= gamma2).sum())
assert np.abs(ct0 - np.rint(ct0)).max() < 1e-6, "FFT lost integer precision"
# Step 4: compute hint probability
prob = np.maximum(0, np.abs(np.rint(ct0)) - beta - 0.5) / denom
n_hints = (rng.random(prob.shape) < prob).sum(axis=(1, 2))
total_hints += int(n_hints.sum())
# Step 5: Number of rejected hints
rejected += int((n_hints > omega).sum())
done += b
est = rejected / done
return total_hints / done, est, math.sqrt(est * (1 - est) / done), ct0_rej
def p_hint_compare(k, tau, d, gamma2, beta, omega, trials, batch=2000, seed=19820116):
"""Monte Carlo comparison of P_hint under the conditioned model (this note)
against the unconditioned model Yan, Luo and Wang use for this term, on
shared (c, t0, u) draws so the comparison is a paired one and not itself
noisy on top of the two individual estimates.
Conditioned: P(hint at i | a) = max(0, a - beta - 1/2) / (2*(gamma2-beta) - 1)
Unconditioned: P(hint at i | a) = min(1, a / (2*gamma2))
Returns (p_cond, p_uncond, se_cond, se_uncond).
"""
rng = np.random.default_rng(seed)
half = 2 ** (d - 1)
denom_cond = 2 * (gamma2 - beta) - 1
denom_uncond = 2 * gamma2
rejected_cond = rejected_uncond = done = 0
while done < trials:
b = min(batch, trials - done)
c = np.zeros((b, 256))
pos = np.argsort(rng.random((b, 256)), axis=1)[:, :tau]
np.put_along_axis(c, pos,
(rng.integers(0, 2, size=(b, tau)) * 2 - 1).astype(float),
axis=1)
t0 = rng.integers(-half + 1, half + 1, size=(b, k, 256)).astype(float)
ct0 = np.real(np.fft.ifft(np.fft.fft(t0 * PSI, axis=2)
* np.fft.fft(c * PSI, axis=1)[:, None, :],
axis=2) * IPSI)
a = np.abs(np.rint(ct0))
prob_cond = np.maximum(0, a - beta - 0.5) / denom_cond
prob_uncond = np.minimum(1, a / denom_uncond)
u = rng.random(a.shape)
n_cond = (u < prob_cond).sum(axis=(1, 2))
n_uncond = (u < prob_uncond).sum(axis=(1, 2))
rejected_cond += int((n_cond > omega).sum())
rejected_uncond += int((n_uncond > omega).sum())
done += b
p_cond = rejected_cond / done
p_uncond = rejected_uncond / done
se_cond = math.sqrt(p_cond * (1 - p_cond) / done)
se_uncond = math.sqrt(p_uncond * (1 - p_uncond) / done)
return p_cond, p_uncond, se_cond, se_uncond
def p_hint_normal_compare(k, tau, d, gamma2, beta, omega, trials, batch=2000, seed=19820116):
"""Compare the true (simulated) P_hint against a normal approximation
that uses N's correct mean and variance (both estimated from the same
trials), to quantify how much the normal approximation underestimates
the tail - the correlation between coefficients of one c*t0 polynomial
(Step 3 in the 'why there is no usable formula' argument) fattens the
true upper tail beyond what a normal with the right variance predicts.
Returns (mean, sigma, p_true, p_normal).
"""
rng = np.random.default_rng(seed)
half = 2 ** (d - 1)
denom = 2 * (gamma2 - beta) - 1
done = 0
sum_n = sum_n2 = 0.0
rejected = 0
while done < trials:
b = min(batch, trials - done)
c = np.zeros((b, 256))
pos = np.argsort(rng.random((b, 256)), axis=1)[:, :tau]
np.put_along_axis(c, pos,
(rng.integers(0, 2, size=(b, tau)) * 2 - 1).astype(float),
axis=1)
t0 = rng.integers(-half + 1, half + 1, size=(b, k, 256)).astype(float)
ct0 = np.real(np.fft.ifft(np.fft.fft(t0 * PSI, axis=2)
* np.fft.fft(c * PSI, axis=1)[:, None, :],
axis=2) * IPSI)
a = np.abs(np.rint(ct0))
prob = np.maximum(0, a - beta - 0.5) / denom
n_hints = (rng.random(prob.shape) < prob).sum(axis=(1, 2))
sum_n += n_hints.sum()
sum_n2 += (n_hints.astype(np.float64) ** 2).sum()
rejected += int((n_hints > omega).sum())
done += b
mean = sum_n / done
sigma = math.sqrt(sum_n2 / done - mean ** 2)
p_true = rejected / done
# continuity-corrected normal tail
z = (omega + 0.5 - mean) / sigma
p_normal = 0.5 * math.erfc(z / math.sqrt(2))
return mean, sigma, p_true, p_normal
def mean_hint_count(k, tau, d, gamma2, beta):
""" Mean hint count E[N] for ML-DSA, analytically.
N is the number of 1 bits in the hint vector h (FIPS 204 Alg 7, line 27).
Its mean follows from linearity of expectation, which holds whether or not
the coefficients are independent -- unlike the tail P(N > omega), which
depends on the joint distribution and has no tractable closed form.
E[N] = sum over coefficients of P(hint at i)
= 256*k * E[ max(0, |c*t0_i| - beta - 1/2) ] / (2*(gamma2-beta) - 1)
The per-coefficient probability is conditioned on line 23 having passed,
which confines r0 to |r0| < gamma2 - beta, and averaged over the sign of
c*t0_i (a fair coin flip given its magnitude, by t0's symmetry). See the
companion note.
A coefficient of c*t0 is a sum of tau terms +-t0_j with t0_j uniform over
2^d values, so it is zero-mean with
sigma = 2^d * sqrt(tau/12) and E|X| = sigma * sqrt(2/pi)
Since beta/sigma is between 0.005 and 0.012 for all parameter sets, |X|
exceeds beta in all but a negligible fraction of cases, so
E[ max(0, |X| - beta - 1/2) ] ~= E|X| - beta - 1/2
"""
sigma = 2 ** d * math.sqrt(tau / 12) # spread of a c*t0 coefficient
e_abs = sigma * math.sqrt(2 / math.pi) # E|X| for a zero-mean Gaussian
per = (e_abs - beta - 0.5) / (2 * (gamma2 - beta) - 1)
return 256 * k * per, sigma, e_abs, per
def run_simulation(trials, params):
"""Run the Monte Carlo for every parameter set and print a summary table.
`params` is the list of MLDSAParams to run (see params.py) - passed in
explicitly rather than defaulted, so this module stays decoupled from
params.py.
"""
phint = dict()
for p in params:
sim, est, se, ct0_rej = p_hint(p.k, p.tau, p.d, p.gamma2, p.beta, p.omega, trials)
print(f"{p.name}: P_hint = {est:.5f} +- {se:.5f} - {sim:.2f} - {ct0_rej}")
phint[p.name] = [sim, est, se]
hdr = (f"{'variant':<12}{'omega':>7}{'sigma':>9}{'E|X|':>9}{'per-coeff':>12}"
f"{'E[N]':>9}{'simulated':>11}")
print(hdr)
print("-" * len(hdr))
for p in params:
lam, sigma, e_abs, per = mean_hint_count(p.k, p.tau, p.d, p.gamma2, p.beta)
print(f"{p.name:<12}{p.omega:7d}{sigma:9.0f}{e_abs:9.0f}{per:12.5f}"
f"{lam:9.2f}{phint[p.name][0]:11.1f}")
return phint
# Running 100M trials per parameter set at site-build time would be
# impolite - it takes on the order of hours. Regenerated after the
# conditioned per-coefficient formula moved from a-beta-1 to a-beta-1/2
# (see @sec-hint) with:
# python3 -m mldsa_rejection.run_phint --trials 20000000 --json result.json
#
# variant trials E[N] sim E[N] analytic P_hint std err seconds
# -------------------------------------------------------------------------------
# ML-DSA-44 20000000 63.07 62.98 0.01458 0.00003 480.1
# ML-DSA-65 20000000 38.23 38.19 0.00380 0.00001 714.7
# ML-DSA-87 20000000 56.75 56.70 0.00785 0.00002 917.6
P_HINT = {
"ML-DSA-44": (63.07, 0.01458, 0.00003),
"ML-DSA-65": (38.23, 0.00380, 0.00001),
"ML-DSA-87": (56.75, 0.00785, 0.00002),
}
df = pd.DataFrame([
{
"**Variant**": p.name,
"**Mean hint count**": f"{P_HINT[p.name][0]:.1f}",
r"$\omega$": p.omega,
r"$P_{\mathrm{hint}}$": f"{P_HINT[p.name][1]:.5f} ± {P_HINT[p.name][2]:.5f}",
}
for p in PARAMS
])
Markdown(df.to_markdown(index=False, disable_numparse=True, colalign=("left", "left", "left", "left")))
