Generate Bessel I0 coefficients in Python and fix cancellations while evaluating the Kaiser-Bessel window. - #248
Draft
jenskeiner wants to merge 2 commits into
Draft
Conversation
![github-code-quality[bot]](https://avatars.githubusercontent.com/u/9919?s=60&v=4){kind=small}
Accuracy report2213 unchanged · 95 improved · 0 regressed
|
jenskeiner
marked this pull request as draft
August 26, 2026 11:00
Contributor
Author
|
Leaving as draft until benchmarks show no regressions. |
jenskeiner
force-pushed
the
feature/bessel
branch
from
August 26, 2026 14:42
fd1f9a0 to
0803eca
Compare
…-Bessel evaluations.
jenskeiner
force-pushed
the
feature/bessel
branch
from
August 26, 2026 14:52
0803eca to
a34d482
Compare
jenskeiner
changed the title
The generated I0 was accurate but slow: 10.8 ns/call against 2.5 for the
hand-written tables it replaced, and phi_hut inherited a factor of six. The
coefficients were never the problem, the evaluation shape was.
bessel_i0.c
- Branch 2 ships as monomial coefficients in u = 1/x rather than a Chebyshev
series in t. A sum can be regrouped; the Clenshaw recurrence cannot, and that
alone cost the asymptotic branch a factor of four.
- Both branches sum as four independent chains over the coefficients sharing a
residue mod 4 (poly4), so the dependent-operation count is a quarter of a
Horner pass. Table lengths are multiples of four, so the loop needs no
prologue.
- exp(x)*sqrt(u) instead of a second divide by sqrt(x), and log I0 = x +
LOG(SQRT(u)*p) instead of two separate logarithms.
Regrouping a sum is only safe when it cannot cancel, and that turned out to pin
the split. P1 is fine at any degree, its coefficients all positive. P2 is not:
its monomial coefficients grow with the degree, and past roughly degree 24 the
growth factor leaves 1 and climbs -- 466 at MANT_DIG 64 with a split of 10,
1.2e9 at 113 with a split of 15. The split now moves per format to the smallest
value holding that factor at 1, and branch 1 absorbs the degree. Long double
splits at 15 and quadruple at 25.
window.c
- phi_hut used the log domain for I0(a)/I0(m b) to keep the peak from
overflowing, which cost a LOG. Adding bessel_i0_exp_scaled -- exp(-x)*I0(x),
in (0, 1] -- lets it be exp(a - m b) * i0e(a) / i0e(m b) instead. Above the
split the exponential cancels against the asymptotic form, so i0e is cheaper
than I0 itself, and a - m b = -m t^2/(ra + b) is still formed without the
cancellation the log domain was there to avoid.
- 1/i0e(m b) is cached per axis, so the call carries no division. Both factors
are bounded, so the peak normalisation cannot overflow at any m b.
Measured (double, aarch64, glibc, ns/call):
bessel_i0 [0,18.85] 2.76 old 6.37 before 2.30 now
[0,94.2] 2.51 old 10.82 before 3.00 now
bessel_i0_log ~4.5 old 15.27 before 4.65 now
kb_phi_hut m=11 4.08 old 23.66 before 5.79 now
kb_phi m=11 3.83 old 4.01 before 3.63 now
Accuracy is unchanged or better. bessel_i0 is 2.31 ulp on (0,100] against 8.20
for the pre-generator tables and 2.39 before this change; kb_phi_hut matches the
log-domain form it replaces to within noise (5.76 vs 5.25 ulp at m=11 sigma=2)
and stays 10-20x clear of the naive I0(a)/I0(m b) ratio.
tests/besselgen
- emulate mirrors poly4, so the fixed-precision emulation still reproduces the
arithmetic the C performs.
- The precision-parametrised tests cover 64 and 113 bits as well, at the
working precision the generator uses. That is the only cover Intel double
extended gets on a host whose long double is binary128.
- New tests for the two invariants the scheme now rests on: P2 cannot cancel at
any pinned split, and every table length is a multiple of four.
Also fixes an unterminated comment in window.c that had swallowed the
kb_inv_peak header.
Claude-Session: https://claude.ai/code/session_01GPi5g9NWB2rgkdr53xN92D
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Replace the Mathematica notebook that produced the I0 coefficient tables with
tools/besselgen, a deterministicmpmathminimax fitter, and rework the evaluation.Generator (
tools/besselgen)kernel/util/bessel_i0_data.hfor all four floating-point formats.I0(kernel/util/bessel_i0.c)x <= splitEvaluate
I0(x) = 1 + y*P1(y),y = (x/2)^2, Horner over positive coefficients. Replaces a Chebyshev-basis rational evaluated at y = x*x, far outside the domain that basis is conditioned for.bessel_i0_logandbessel_i0_logtailuseLOG1P, i.e.log(I0(x)) = LOG1P(y*P1(y))which avoids cancellations wherelog I0(x) ~ x^2/4near zero.x > splitEvaluate
I0(x) = exp(x) / sqrt(x) * p. Forpmap(split, inf)to(-1,1)and use Clenshaw. Use a slightly different order of evaluation ofexp(x)alone would already overflow. This correctly evaluates results whereI0(x)itself is still representable.Kaiser-Bessel window (
kernel/util/window.c)Phi_hut now always takes the log form. Several other improvements, e.g. use
(b-t)(b+t)and(m-nx)(m+nx)instead of differences of squares, andEXPM1instead of a difference of exponentials for sinh.Tests
kaiser_bessel_cancellation case: phi_hut against 60-digit reference values for the near-peak bins and the band edge.bessel_i0bounds tightened from the measured worst case: 24 -> 8 single, 58 -> 16 quadruple.See bessel-i0-kaiser-bessel-accuracy.html for more details.