TST: kelvin: add a Wronskian identity test - #236
Open
EazyHood wants to merge 3 commits into
Open
Conversation
x * (Be*Kep - Bep*Ke) == -1, with Be = ber + i*bei, Ke = ker + i*kei and Bep, Kep their derivatives -- exactly what xsf::kelvin() returns. One check covers all eight functions, and it needs no reference table, so it stays valid whatever is decided about the |x| < 10 branch point in scipygh-232. Verified against mpmath at 50 digits over x in [0.05, 25]: largest deviation 1.4e-51. On the current implementation the worst residual over the 20 sampled points is 8.2e-11, at x = 9.9.
github-actions
Bot
added
the
Maintenance
The list fits on one 117-column line; splitting it at 70 was below the project's 120-column limit and clang-format packs it.
j-bowhay
reviewed
Aug 4, 2026
Comment on lines
+15
to
+16
| // Verified against mpmath at 50 digits for x in [0.05, 25]: the largest deviation | ||
| // is 1.4e-51, i.e. mpmath's own rounding. |
Member
There was a problem hiding this comment.
What does this comment refer to?
j-bowhay
reviewed
Aug 4, 2026
Comment on lines
+17
to
+20
| // | ||
| // This needs no reference table -- the eight functions are checked against each | ||
| // other -- so it stays valid whatever is decided about the |x| < 10 branch point | ||
| // discussed in gh-232. |
Member
There was a problem hiding this comment.
Suggested change
| // | |
| // This needs no reference table -- the eight functions are checked against each | |
| // other -- so it stays valid whatever is decided about the |x| < 10 branch point | |
| // discussed in gh-232. |
The mpmath line described a separate check, not anything this test does, and the scipygh-232 note was removed per the review suggestion.
Author
|
Both removed. The mpmath line described a separate check I ran while writing the test, not The workflows on this PR are all sitting in action_required and have never run, so there is no |
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
2 participants
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.
Follows up on the Wronskian test offered in gh-232, which @steppi said he'd be happy to see.
What it checks
The Kelvin functions satisfy
with
Be = ber + i*bei,Ke = ker + i*kei, andBep,Keptheir derivatives — exactly the four valuesxsf::kelvin()returns. Its real part is(ber*kerp - berp*ker) - (bei*keip - beip*kei)and its imaginary part is(ber*keip - berp*kei) + (bei*kerp - beip*ker), so a single check exercises all eight functions.I derived it by searching bilinear combinations against
mpmathrather than quoting it: overxin[0.05, 25]at 50 digits the largest deviation is 1.4e-51, which is mpmath's own rounding.It needs no reference table — the functions are checked against each other. That is why it does not depend on how the
|x| < 10branch point question in gh-232 is resolved.What it measures on the current implementation
Sweeping 52 points against
scipy1.17.1:The worst point in the whole sweep sits just below the
|x| < 10branch point, and crossing it improves by ~111x. That is the claim gh-232 makes, arrived at without mpmath and without any reference values of mine — it is the library disagreeing with itself.On the tolerance
atol = 1e-9passes today with about two orders of margin, so this lands as a non-regression test rather than a failing assertion. At1e-12it fails atx = 9.9on currentmain. I chose the looser bound deliberately so the test is mergeable independently of the branch-point work; if you would rather it fail loudly and block until that is fixed, say so and I will tighten it.What I have not done
I have not compiled this. Catch2 3, Arrow and Parquet are not available on the machine I measured from, so the numbers above are verified but the build is not. It goes in
tests/xsf_tests/rather thantests/scipy_special_tests/because it needs no parquet tables, andCMakeLists.txtpicks it up through the existingfile(GLOB TEST_SOURCES "*/test_*.cpp")— no build-system change. If CI finds a compile error I will fix it promptly; the one thing I would guess at is whetherCAPTURE()printsstd::complex<double>, which is diagnostic only.AI disclosure: prepared with AI assistance (Claude Opus 5 via Claude Code). The identity and every figure above were measured, and I am responsible for the change as submitted.