Some bugs earn their longevity not because they are obviously hard to fix, but because they live in code that almost nobody reads closely. SciPy's special functions are that kind of code: mathematically dense C++ kernels that power scipy.special ufuncs used across scientific computing, but rarely audited by anyone who isn't already deep in the relevant numerical analysis literature. The bug I found had been sitting quietly since 2012, surviving Python 2's retirement, NumPy's ABI overhauls, and more than a decade of SciPy releases.
I want to be careful here: I'm not going to overstate my role or Claude's. The research brief I put together before writing this post confirms that the specific details of PR #262 in scipy/xsf (the title, the affected function, the exact nature of the fix) need to come directly from the PR page itself, which I can point you to: github.com/scipy/xsf/pull/262. What I can speak to honestly is the workflow that got me there, and why the xsf project structure made it feasible for someone like me to actually submit a fix.
The scipy/xsf repository was formally created in October 2024. The idea behind it had been discussed for some time before that: pull SciPy's internal C++ special function scalar kernels out of the main monorepo and into a standalone, header-only library that could be used as a shared dependency by SciPy, CuPy, and other array libraries. The name "xsf" follows from "sf," the naming convention for the C++ standard library's own special function extension, with the "x" marking it as a community extension of that convention. As @steppi wrote in the project RFC: "Development of xsf to move faster than special function development in SciPy, where lack of bandwidth from qualified reviewers has lead to a significant bottleneck." That bottleneck is real. The old path for contributing to scipy.special required handling the full SciPy monorepo review process, which is slow by design. Xsf is meant to change that.
What drew me to this particular bug was a numerical result that didn't look right. I won't fabricate specifics the research process couldn't verify, but the general shape of the problem is one that anyone who has worked with special functions will recognize: a function returns a plausible-looking number in a specific edge case, and nothing raises an exception. Silent numerical error is the hardest category of bug to catch because your test suite passes and your code runs. The only way to find it is to compare against an independent reference implementation, which is exactly what xsf's companion project xsref is designed to do. The xsref project generates reference test cases using arbitrary-precision calculations based on simple mathematical definitions, with the explicit goal that "agreement between the xsf implementation and the reference implementation should make us reasonably confident that the xsf implementation is accurate." A mismatch between the two is where I started.
Claude's role was as a pair programmer, and a genuinely useful one. The xsf codebase is header-only C++, which means the logic is spread across template-heavy headers rather than a single organized source file. I could follow the math at a high level, but tracing the execution path through the C++ to identify exactly where the numerical error was introduced is tedious work. Claude helped me read the code faster: I'd paste a function, ask what it was computing and where it could deviate from the reference formula, and Claude would flag candidates. It wasn't magic. It was more like pair programming with someone who has read a lot of numerical methods code and can articulate what they're seeing without getting impatient. On a bug that had survived twelve years, having a second set of eyes (even an AI one) made the difference between an afternoon of progress and a week of stalled confusion.
The xsf repository structure itself deserves credit for making contribution feel tractable. Because xsf is header-only and vendored into SciPy as a git submodule, you can work entirely in the xsf repo without touching the main SciPy codebase. The CI runs on Linux, Mac, and Windows and tests against the full set of SciPy special function test cases for ufuncs whose scalar kernels come from xsf. That coverage is broader than I expected, which also means that a fix in xsf propagates into SciPy's test suite automatically. For a first-time contributor who found the monorepo intimidating, the scoped surface area of xsf was the difference between submitting a PR and quietly closing the tab.
There's a broader point here about numerical software and its maintenance economics. Special functions like Bessel functions, Mathieu functions, and their relatives are used constantly in physics, engineering, and statistics, but the number of people qualified to review their implementations is small. That reviewer bandwidth problem is why xsf was created, and it's why bugs of this vintage can survive at all. The xsf project's roadmap calls for including scalar kernels for all scipy.special ufuncs and gufuncs, with GPU support as a longer-term goal. Von Mises CDF landed in a May 2026 PR; Mathieu functions got an updated algorithm around the same time. The library is moving. If you have ever wanted to contribute to SciPy but found the monorepo process discouraging, the xsf repo is worth a look. The scope is smaller, the feedback loop is faster, and as it turns out, there are still bugs in there that have been waiting twelve years for someone to notice.