Skip to content

Fix Jacobi-Anger degree search stopping at a Bessel zero crossing #1911

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Sanjays2402 wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Fix Jacobi-Anger degree search stopping at a Bessel zero crossing #1911

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 12 additions & 4 deletions qualtran/linalg/polynomial/jacobi_anger_approximations.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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))))

@gemini-code-assist gemini-code-assist Bot Jul 25, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

high

The minimum degree d_min is calculated using float(t). However, if t is negative, float(t) will be negative, resulting in d_min = 1. Since |J_n(t)| = |J_n(|t|)|, the Bessel function magnitude is symmetric with respect to t, and the oscillatory region still extends up to |t|. For negative t, this would cause the search to start at 1 and potentially truncate at a zero-crossing, re-introducing the bug. Using abs(float(t)) ensures that the decaying region is correctly identified for both positive and negative values of t.

Suggested change
d_min = max(1, int(np.ceil(float(t))))
d_min = max(1, int(np.ceil(abs(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


Expand Down
18 changes: 18 additions & 0 deletions qualtran/linalg/polynomial/jacobi_anger_approximations_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,24 @@ def test_exp_cos_approximation(t: float, precision: float):
)


@pytest.mark.parametrize("t", [20, 50])

@gemini-code-assist gemini-code-assist Bot Jul 25, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

high

To ensure that negative values of t are also correctly handled and do not suffer from the same zero-crossing truncation issue, we should include negative test cases in the parameter list.

Suggested change
@pytest.mark.parametrize("t", [20, 50])
@pytest.mark.parametrize("t", [-50, -20, 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):
Expand Down