Version: qualtran 0.7.0 (pip), Python 3.14, scipy 1.18 / numpy 2.x, macOS. File: qualtran/linalg/polynomial/jacobi_anger_approximations.py.
Minimal repro
from qualtran.linalg.polynomial.jacobi_anger_approximations import (
degree_jacobi_anger_approximation,
)
degree_jacobi_anger_approximation(150.0, precision=5e-4) # returns 1 (should be ≈160+)
degree_jacobi_anger_approximation(1500.0, precision=5e-4) # returns 1529 (correct)
degree_jacobi_anger_approximation(225000.0, precision=5e-4) # returns 1 (should be ≈225,000+)
Some inputs instead crash on the function’s own final assert (various τ at precision=1e-3).
How often
Self-contained scan — paste and run; reproduces the table below exactly (verified on a clean qualtran==0.7.0 install, 2026-08-16):
import numpy as np
from qualtran.linalg.polynomial.jacobi_anger_approximations import (
degree_jacobi_anger_approximation as deg)
taus = 10 ** np.random.default_rng(20260810).uniform(1, 6, 4000)
for eps in (1e-3, 5e-4, 1e-4):
wrong = crash = 0
for t in taus:
try:
# |J_n(t)| only decays for n > t, so any d < t/2 cannot be sound
if deg(float(t), precision=eps) < t / 2:
wrong += 1
except AssertionError:
crash += 1
print(f"eps={eps:<7} silently wrong {100wrong/len(taus):5.1f}% "
f"AssertionError {100crash/len(taus):4.1f}%")
Log-uniform scan of 4,000 values of t ∈ [10, 10⁶]:
| precision |
silently wrong (d < t/2) |
AssertionError |
median underestimate factor |
| 1e-3 |
40.4% |
6.4% |
~2800× |
| 5e-4 |
35.5% |
0% |
~630× |
| 1e-4 |
13.5% |
0% |
~170× |
(The d < t/2 flag is deliberately conservative — it only counts degrees that cannot possibly be sound. Rates are stable across seeds: an independent seed at 2,000 samples gives 38.4% / 33.7% / 12.8%.)
Root cause
The truncation criterion |J_{d+1}(t)| ≤ ε is only meaningful in the decay regime n > t. For n ≲ t, J_n(t) oscillates through zeros with envelope ≈ √(2/πt), and the implementation makes three assumptions that fail there:
1. The doubling probe tests a single value. term_too_small(d) checks one |J_d(t)|. If the doubling sequence 1, 2, 4, … lands near an oscillation zero, the loop exits with d far inside the oscillatory band, where thousands of O(t^{-1/2}) terms remain in the tail.
2. bisect.bisect_left(range(d), True, key=term_too_small) requires the predicate to be monotone (False…False,True…True). In the oscillatory band it is not, so bisect returns an arbitrary crossing; the final assert not term_too_small(d) and term_too_small(d+1) only checks a local crossing, which any oscillation zero satisfies — so it usually passes on wrong answers, and occasionally fails outright (the AssertionError mode).
3. Deterministic total failure at large t: once √(2/πt) < ε — i.e. t > 2/(πε²), about 6.4×10⁵ at ε = 1e-3 — every n in the oscillatory band tests “too small,” so the first probe succeeds and the function returns ~1 for essentially all such t.
Impact
HamiltonianSimulationByGQSP consumes this degree with τ = t·α (block-encoding 1-norm times time). τ in the 10⁴–10⁶ range is routine for chemistry/materials-scale estimates, so resource counts derived through this path can silently underestimate T/Toffoli cost by 2–3 orders of magnitude. (Found while cross-validating qualtran against other resource estimators; at τ ≈ 150 on a 6×6 XXZ evolution instance the GQSP cost came out ~1000× low. It has since fired in a second, unrelated realistic costing: a 2D transverse-field Ising dynamics instance at τ = 960, precision 5e-4, returned degree 31 — a ~32× degree underestimate.)
Suggested fix
Restrict the search to the monotone decay regime: start both the doubling probe and the bisect at n₀ = ceil(|t|) (|J_n(t)| is monotonically and super-exponentially decreasing in n for n > |t|), e.g.
n0 = max(1, math.ceil(abs(t)))
d = n0
while not term_too_small(d):
d *= 2
d = n0 + bisect.bisect_left(range(n0, d), True, key=term_too_small) - 1
This also fixes the AssertionError mode (the predicate is monotone on [n₀, ∞)) and matches the symbolic branch’s asymptotic d = O(t + log(1/ε)/log log(1/ε)). A stricter follow-up would bound the tail sum rather than a single term, but the n₀ floor is the load-bearing correction. Happy to submit a PR with the fix + a regression test over the scanned τ range if useful.
Version: qualtran 0.7.0 (pip), Python 3.14, scipy 1.18 / numpy 2.x, macOS. File: qualtran/linalg/polynomial/jacobi_anger_approximations.py.
Minimal repro
from qualtran.linalg.polynomial.jacobi_anger_approximations import (
degree_jacobi_anger_approximation,
)
degree_jacobi_anger_approximation(150.0, precision=5e-4) # returns 1 (should be ≈160+)
degree_jacobi_anger_approximation(1500.0, precision=5e-4) # returns 1529 (correct)
degree_jacobi_anger_approximation(225000.0, precision=5e-4) # returns 1 (should be ≈225,000+)
Some inputs instead crash on the function’s own final assert (various τ at precision=1e-3).
How often
Self-contained scan — paste and run; reproduces the table below exactly (verified on a clean qualtran==0.7.0 install, 2026-08-16):
import numpy as np
from qualtran.linalg.polynomial.jacobi_anger_approximations import (
degree_jacobi_anger_approximation as deg)
taus = 10 ** np.random.default_rng(20260810).uniform(1, 6, 4000)
for eps in (1e-3, 5e-4, 1e-4):
wrong = crash = 0
for t in taus:
try:
# |J_n(t)| only decays for n > t, so any d < t/2 cannot be sound
if deg(float(t), precision=eps) < t / 2:
wrong += 1
except AssertionError:
crash += 1
print(f"eps={eps:<7} silently wrong {100wrong/len(taus):5.1f}% "
f"AssertionError {100crash/len(taus):4.1f}%")
Log-uniform scan of 4,000 values of t ∈ [10, 10⁶]:
(The d < t/2 flag is deliberately conservative — it only counts degrees that cannot possibly be sound. Rates are stable across seeds: an independent seed at 2,000 samples gives 38.4% / 33.7% / 12.8%.)
Root cause
The truncation criterion |J_{d+1}(t)| ≤ ε is only meaningful in the decay regime n > t. For n ≲ t, J_n(t) oscillates through zeros with envelope ≈ √(2/πt), and the implementation makes three assumptions that fail there:
1. The doubling probe tests a single value. term_too_small(d) checks one |J_d(t)|. If the doubling sequence 1, 2, 4, … lands near an oscillation zero, the loop exits with d far inside the oscillatory band, where thousands of O(t^{-1/2}) terms remain in the tail.
2. bisect.bisect_left(range(d), True, key=term_too_small) requires the predicate to be monotone (False…False,True…True). In the oscillatory band it is not, so bisect returns an arbitrary crossing; the final assert not term_too_small(d) and term_too_small(d+1) only checks a local crossing, which any oscillation zero satisfies — so it usually passes on wrong answers, and occasionally fails outright (the AssertionError mode).
3. Deterministic total failure at large t: once √(2/πt) < ε — i.e. t > 2/(πε²), about 6.4×10⁵ at ε = 1e-3 — every n in the oscillatory band tests “too small,” so the first probe succeeds and the function returns ~1 for essentially all such t.
Impact
HamiltonianSimulationByGQSP consumes this degree with τ = t·α (block-encoding 1-norm times time). τ in the 10⁴–10⁶ range is routine for chemistry/materials-scale estimates, so resource counts derived through this path can silently underestimate T/Toffoli cost by 2–3 orders of magnitude. (Found while cross-validating qualtran against other resource estimators; at τ ≈ 150 on a 6×6 XXZ evolution instance the GQSP cost came out ~1000× low. It has since fired in a second, unrelated realistic costing: a 2D transverse-field Ising dynamics instance at τ = 960, precision 5e-4, returned degree 31 — a ~32× degree underestimate.)
Suggested fix
Restrict the search to the monotone decay regime: start both the doubling probe and the bisect at n₀ = ceil(|t|) (|J_n(t)| is monotonically and super-exponentially decreasing in n for n > |t|), e.g.
n0 = max(1, math.ceil(abs(t)))
d = n0
while not term_too_small(d):
d *= 2
d = n0 + bisect.bisect_left(range(n0, d), True, key=term_too_small) - 1
This also fixes the AssertionError mode (the predicate is monotone on [n₀, ∞)) and matches the symbolic branch’s asymptotic d = O(t + log(1/ε)/log log(1/ε)). A stricter follow-up would bound the tail sum rather than a single term, but the n₀ floor is the load-bearing correction. Happy to submit a PR with the fix + a regression test over the scanned τ range if useful.