If I am not mistaken, which is extremely easy with this file, these are all the public functions we use. I tested half of them by Wolfram Alpha by reading what it does from the docstring, but some I don't know much about the subject and probably we'll rely on the existing tests.

All tests are now passing on my local machine. However, of course because C, occasionally I'm getting a windows access violation crash. I don't know which one is the offending one but I'll keep digging. If anybody knows how to pinpoint the actual test that originates the violation that would be awesome (maybe a tool I don't know yet exists in Linux or Mac?).

Other than that, I tested most of the if branches I see in the functions and when succeeds it passes all so let's see what CI would say.

Every commit is one specific file hence we can keep the commit structure just in case.

Nevermind CI actually caught them correctly. I'll fix those.

I don't understand why cdflib test is failing on this job, we just merged the same code with all passing.

I don't understand why cdflib test is failing on this job, we just merged the same code with all passing.

Is it possible the CI for #8404 last ran without properly merging main, and thus didn’t have those tests?

I've confirmed locally that that the test fails with tol = 1e-15 but not with tol = 1e-8. I think #8404 should be reverted.

I think I understand what's going on, and why the tolerance was chosen at the level it was. In the following, I will refer to the current estimate in a root finding iterative process as $x_{\text{current}}$ and the previous estimate as $x_{\text{previous}}$.

For root finding algorithms with better than linear convergence, $|x_{\text{current}} - x_{\text{previous}}|$ will lag behind the actual error, because the actual root should be closer to $x_{\text{current}}$ than to $x_{\text{previous}}$. If it wasn't closer, you couldn't possibly have superlinear convergence. For algorithm R from Bus and Dekker 1975, used in cdflib, $x_{\text{current}}$ and $x_{\text{previous}}$ bracket the root, so

$$|x_{\text{current}} - x_{\text{previous}}| = |x_{\text{current}} - \text{root}| + |\text{root} - x_{\text{previous}}|$$

This means that the value being compared against tol at each iteration is not an estimate of the current error, it is the sum of the current error and the previous error, which for fast enough convergence will be an estimate of the previous error.

The rate of convergence of algorithm R used in cdflib is bounded below by 1.6202 (but could be faster depending on properties of the function). This made me think that

10**(np.log10(np.finfo(float).eps)/1.6202) = 2.1802262222980188e-10

might be a better tolerance. This takes into account the rate of convergence to estimate what $|x_{\text{current}} - x_{\text{previous}}|$ would be if $|x_{\text{current}} - \text{root}| = \epsilon$

I've found the test still fails in that case. I think because 1.6202 is only a lower bound on the rate of convergence, thus it's still possible to make additional steps after optimal convergence has been reached, causing accumulation of rounding errors to swamp the benefit of additional iterations.

@ilayn, while your at it, can you fix #17681? See #11740 for the functions that were removed, and cross reference against what's supposed to be in the public API according to the docs.

I checked and the only wrapper missing in _specfun.pyx now is the one for pbvv. We should probably add a test for pbvv_seq,

FWIW, I just did a coverage run on the main branch via python dev.py test --coverage --submodule special (just covers Python code from what I see) and went manually through the _basic.py file which calls a lot of specfun functions. The following functions are not covered:

• clqmn
• clpn
• lamv
• clqn
• pbvv

pbvv is already on your radar from what I see but for the others I am not sure if the newly added tests would cover them. Feel free to ignore if I missed something in the long discussion above.

We should probably add a test for pbvv_seq

If you have known values for the test, I can add it right away.

We should probably add a test for pbvv_seq

If you have known values for the test, I can add it right away.

Interesting, pbvv is almost untested because the mpmath function that's supposed to be the same is different. I'd need to figure out how pbvv relates to the function pcfv in mpmath.

Ah mpmath uses $a$ and pbvv uses $v = -(a + \frac{1}{2})$

Ok, if I haven't made any mistakes, pbvv_seq(2, 3) should approximately equal

(
    np.array([2.976319645712036, 1.358840996329579, 0.5501016716383508]),
    np.array([3.105638472238475, 0.9380581512176672, 0.533688488872053]),
)

we should also fix test_pcfv in test_mpmath.py by using the relation v = -(a + 0.5)

Thanks @dschmitz89, I'll come up with test cases for those others too.

I've just seen that mpmath pbvv test still fails in 27 out of 746 some cases even after adjusting the parameterization (and setting rtol to a generous 1e-3). It may be due to issues in the original implementation or perhaps because mpmath uses a different branch in some parts of the complex plane. It seems lqmn is also untested in the complex plane because mpmath uses a different branch for some points.

For clpn we just need to test that lpn(n, z) is equal to clpmn(0, n, z) for some complex z, since the former is a special case of the later. Can just pick some random test case.

@steppi I first repeated all your comments but the rswfo -10 to 10 correction revealed a segfault on windows that I can't figure out yet because everything started to look the same. After fixing that we can add all the tests.

Thank you very much for your thorough review.

@steppi I first repeated all your comments but the rswfo -10 to 10 correction revealed a segfault on windows that I can't figure out yet because everything started to look the same. After fixing that we can add all the tests.

Thank you very much for your thorough review.

What input is causing the segfault? And happy to review. It's awesome to see this getting done.

The two tests test_obl_rad2 and test_obl_rad2_cv (in test_basic.py) are the two tests lead to segfault on windows. Everything else passes.

The two tests test_obl_rad2 and test_obl_rad2_cv (in test_basic.py) are the two tests lead to segfault on windows. Everything else passes.

Thanks. Also this is really confusing:

A holdover from when scipy.special was just cephes in the very earliest days.

@steppi if we are happy with its state should we merge and do a complete sweep on testing including cdflib ones too on a separate PR? I'm a bit worried that there is too much happening in this one and we probably need a bit more focus on tests alone as @dschmitz89 and @inkydragon noticed already.

@steppi if we are happy with its state should we merge and do a complete sweep on testing including cdflib ones too on a separate PR? I'm a bit worried that there is too much happening in this one and we probably need a bit more focus on tests alone as @dschmitz89 and @inkydragon noticed already.

Yeah, I think this is in good shape. I'm going to go everything quickly one more time and plan to merge today if I find no further issues.

Thanks all for the input above and beyond! Very happy that this is finally happening. Let's finish it off by completing the testing suite.

Quite the milestone here.