From eb331fae7cef65fafeb34b1285c6907a17b9f0fd Mon Sep 17 00:00:00 2001
From: Sanjay Santhanam <51058514+Sanjays2402@users.noreply.github.com>
Date: Sat, 25 Jul 2026 08:12:27 -0700
Subject: [PATCH] Fix Jacobi-Anger degree search stopping at a Bessel zero
 crossing

degree_jacobi_anger_approximation searched for the smallest n with
|J_n(t)| <= precision starting from n = 1. For n < t, |J_n(t)| oscillates
and crosses zero repeatedly, so with a loose precision the doubling loop
and the following bisection could land on a zero crossing rather than the
asymptotic tail: for t = 50, precision = 1e-2 it returned degree 15
(|J_16(50)| = 0.005) instead of 56, and the resulting truncated expansion
was not an approximation of exp(i t cos(theta)) at all.

|J_n(t)| decays monotonically only for n >= t, so the search now starts at
ceil(t) and treats every smaller n as not-yet-small, which keeps the
bisection predicate monotonic.

Fixes #1848
---
 .../polynomial/jacobi_anger_approximations.py  | 16 ++++++++++++----
 .../jacobi_anger_approximations_test.py        | 18 ++++++++++++++++++
 2 files changed, 30 insertions(+), 4 deletions(-)

diff --git a/qualtran/linalg/polynomial/jacobi_anger_approximations.py b/qualtran/linalg/polynomial/jacobi_anger_approximations.py
index 07a84e70c8..c9461f6d31 100644
--- a/qualtran/linalg/polynomial/jacobi_anger_approximations.py
+++ b/qualtran/linalg/polynomial/jacobi_anger_approximations.py
@@ -56,13 +56,21 @@ def degree_jacobi_anger_approximation(t: SymbolicFloat, *, precision: SymbolicFl
     def term_too_small(n: int) -> bool:
         return bool(np.isclose(scipy.special.jv(n, t), 0, atol=float(precision)))
 
-    d = 1
-    while not term_too_small(d):
+    # `|J_n(t)|` oscillates and crosses zero many times for `n < t`, and only decays
+    # monotonically for `n >= t`. Restricting the search to the decaying region keeps
+    # `term_too_small` monotonic, so we do not truncate at a zero crossing.
+    d_min = max(1, int(np.ceil(float(t))))
+
+    def tail_term_too_small(n: int) -> bool:
+        return n >= d_min and term_too_small(n)
+
+    d = d_min
+    while not tail_term_too_small(d):
         d *= 2
 
     # find the smallest `n` such that J_n(z) is too small
-    d = bisect.bisect_left(range(d), True, key=term_too_small) - 1
-    assert not term_too_small(d) and term_too_small(d + 1)
+    d = bisect.bisect_left(range(d), True, key=tail_term_too_small) - 1
+    assert term_too_small(d + 1)
     return d
 
 
diff --git a/qualtran/linalg/polynomial/jacobi_anger_approximations_test.py b/qualtran/linalg/polynomial/jacobi_anger_approximations_test.py
index 166e094140..edf273e1f4 100644
--- a/qualtran/linalg/polynomial/jacobi_anger_approximations_test.py
+++ b/qualtran/linalg/polynomial/jacobi_anger_approximations_test.py
@@ -37,6 +37,24 @@ def test_exp_cos_approximation(t: float, precision: float):
     )
 
 
+@pytest.mark.parametrize("t", [20, 50])
+@pytest.mark.parametrize("precision", [1e-2, 1e-3])
+def test_exp_cos_approximation_loose_precision(t: float, precision: float):
+    """`|J_n(t)|` oscillates through zero for `n < t`, so the degree search must not stop there."""
+    random_state = np.random.RandomState(42 + int(t))
+
+    degree = degree_jacobi_anger_approximation(t, precision=precision)
+    P = np.polynomial.Polynomial(approx_exp_cos_by_jacobi_anger(t, degree=degree))
+    theta = 2 * np.pi * random_state.random(1000)
+    e_itheta = np.exp(1j * theta)
+    np.testing.assert_allclose(
+        P(e_itheta) * e_itheta ** (-degree),
+        np.exp(1j * t * np.cos(theta)),
+        atol=precision * 10,
+        rtol=0,
+    )
+
+
 @pytest.mark.parametrize("t", [2, 3, 5, 10])
 @pytest.mark.parametrize("precision", [1e-5, 1e-7, 1e-10])
 def test_exp_sin_approximation(t: float, precision: float):