ENH: Updated and refactored Mathieu function implementation (part 1).

[steppi]

Reference issue

Partially supersedes #99

What does this implement/fix?

This PR is based on @brorson's PR #99, which gives improved implementations of the Mathieu functions. Algorithmically, #99 is excellent work. This is a promised refactor in order to put the work in #99 into the form that would be most useful for xsf and SciPy. The primary departure is that this PR does not give full implementations of Mathieu functions; it only contains the core numerical kernels that would be useful in either a CPU implementation for SciPy or a GPU implementation for CuPy. The primary functions are for generating the recurrence matrices and for summing the Fourier series. Generating the Fourier coefficients requires finding eigenvectors of symmetric tridiagonal recurrence matrices. #99 used LAPACK's dstevd for this, which is out of scope for xsf. I do not rule out eventually writing a CPU+GPU capable symmetric tridiagonal eigenvector solver, but even then, it would still not be particularly useful for SciPy and CuPy to have full scalar implementations in xsf. For SciPy implementations, it is best to use the stateful-functor with in-ufunc caching pattern established here scipy/scipy#24894, in order to allow reuse of the same computed fourier coefficients across varying values of x. For CuPy, we will want to use the gufunc expansion + reduction pattern proposed here cupy/cupy#9839.

The vision I'm coming to is that xsf will contain simple numerical kernels which will work on both CPU and GPU. If a function is simple enough for the same scalar implementation to work on both, then a full implementation belongs in xsf. Otherwise, xsf should contain the kernels that can be common to both platforms and CPU specific implementations can live in SciPy and GPU specific implementations in CuPy.

This PR only contains kernels used for the characteristic value functions mathieu_a and mathieu_b and the Mathieu functions mathieu_cem and mathieu_sem from SciPy. I plan to address mathieu_modcem1, mathieu_modcem2, mathieu_modsem1, and mathieu_modsem2 in an eventual follow-up. This is primarily to cut down the scope of work in this PR. There are also some redundant computations in the Fourier-Bessel summations from #99 that we can eliminate by pushing through #50, which is worth doing on its own as well.

I have used extensive templating to cut down on the duplicate code in #99. The implementations here should otherwise be completely faithful to #99, with the exception that I have added an AngleUnits policy template arg to the Fourier series summation so that the angular units can be taken to be either degrees or radians. For backwards compatibility, the companion PR to SciPy still takes arguments in degrees. There is an easy path forward to adding a kwarg to the SciPy functions to select between radians and degrees.

The tests here test only the recurrence matrix generation and the fourier summation. The tests for fourier summation use coefficients generated by the SciPy implementation and compare the evaluated sums to reference implementations of Mathieu functions from Mathematica. I've discovered the Wolfram Engine Free Developers program, and used the Python bindings to get reference values. This seems like something that will be useful for xsref too, which has been in a languishing state, but which I still have high hopes for.

I've removed the old Mathieu tests, which were self tests against the inaccurate SciPy code anyway.

Additional information

After cleaning up the stuff adding BLAS and LAPACK to the dependencies, the rebuild of pixi.toml failed. This is because the cupy-tests environment was included in solve-group = "default", but because the CuPy tests are Linux only and pull down things not available on other platforms, this caused solver errors. I simply removed the CuPy environment from that solve group.

AI Generation Disclosure

The test code in xsf/tests/xsf_tests/test_mathieu.cpp was largely AI generated with Google Gemini 3.1 Pro. I checked carefully that everything is correct. The implementation code is all good-old-fashioned artisanal hand-crafted code; this was too fun to delegate to AI.

[steppi]

Oh wow, GPU CI actually ran after I updated the ubuntu runner to 24.04. Tests are failing, but I think will work if we update CuPy in the pixi.lock file. I'll make a separate PR for that.

[steppi]

The formatting errors seem to be because the clang-format version was updated when I updated the pixi.lock file. Let's leave them here and fix them in another PR.

[brorson]

Great to see this get merged in!

[fbourgey]

@steppi, this a first pass.

[fbourgey]

If we think of the very negative $q$ case, should we not use if (std::abs(q) > 1.0) instead?

[steppi]

I think you may be right. @brorson, can you comment on this?

[brorson]

I think that's right. The point here is that we need more terms in the sum when abs(q) > 1. So, yes, if (std::abs(q) > 1.0) is correct. Don't forget to use std::log10(std::abs(q)) on line 167 too.

[fbourgey]

I understand that if X is an eigenvector, then so is -X but can you detail the rule "if the coefficients mostly sum negative, we flip the results" implemented here? I guess the coefficients of X(l) can have any sign in this framework (unless I missed something), so if the sum is very close to 0, can this be unstable? is there a better way?

[steppi]

This is a very good question and I think you've unearthed a potential bug. dstevd will return the eigenvectors normalized to have $l^2$ norm equal to 1, but the chosen sign is arbitrary, so we could either get the eigenvector $\mathbf{x}$ or $-\mathbf{x}$. These functions satisfy normalization conditions $\mathrm{cem}(m, q, 0) > 0$ and $\mathrm{sem}^{'}(m, q, 0) > 0$ (see details section below). $\mathrm{cem}(m, q, 0)$ is equal to the sum of the coefficients, so @brorson's idea here is to use this fact to get the correct sign for $\mathbf{x}$. But the correct normalization condition is $\mathrm{sem}^{'}(m, q, 0) > 0$ with $\mathrm{sem}(m, q, 0)$ identically 0. Rather than the sum of fourier coefficients, the value of this derivative at $x=0$ is

$$\mathrm{sem}^{'}(m, q, 0) = \frac{\mathrm{d}}{\mathrm{d}x}\left[\sum_{k}B_k \sin(kx)\right]_{x=0} = \sum_{k}kB_k$$

so I think to handle FuncParity odd, the above should be replaced with:

        if constexpr (FuncParity == Parity::Even) {
                x += X(l);
            } else {
                double r = sqrt_di<Parity::Odd, OrderParity>(l);
                x += r * X(l);
            }

@brorson, does that sound reasonable?

Details

For $q = 0$ we have (https://dlmf.nist.gov/28.2#E29):

$$\mathrm{cem}(m, 0, x) = \cos(mx)$$
$$\mathrm{sem}(m, 0, x) = \sin(mx)$$

which implies that the sign conventions I mentioned above hold for $q = 0$. DLMF glosses over the convention for general $q$ by saying:

the ambiguity of sign being resolved by (28.2.29) when q=0 and by continuity for the other values of q.

It wasn't obvious to me at first, but suppose for instance there is some $q$ such that $\mathrm{cem}(m, q, 0) < 0$, by continuity, there must be some $q^{\star}$ for which $\mathrm{cem}(m, q^{\star}, 0) = 0$. Since we know $\mathrm{cem}(m, q^{\star}, x)$ is an even function, it follows that $\mathrm{cem}^{'}(m, q^{\star}, 0) = 0$. The only solution to Mathieu's equation satisfying these boundary conditions must be identically 0, from which we could derive a contradiction.

[steppi]

As for stability, the sum of coefficients is 0 at roots, and calculating these Mathieu functions seems to be inherently unstable near roots (or roots of the derivative). I found even the arbitrary precision Mathematica reference struggles near roots. We still have a major improvement over Zhang and Jin, so I think it's fine if there's still some rough edges.

The above was nonsense. Sorry. The sum of coefficients is zero implies mathieu_cem is identically zero for all x since the function and it's derivative would be 0 at x = 0 and $f(x) = 0$ is the only solution of Mathieu's equations satisfying these boundary conditions.

I don't think this situation can occur for any values of m or q.

[brorson]

Yes, I think sum(k*B[k]) is reasonable. Go ahead and change it. I hope to have students run my test suite on the results at some point so if your idea doesn't work then we will pick it up.

[steppi]

Thanks @fbourgey for the excellent review so far. I will fix all of the typos. I've tried to clarify the eigenvector sign situation, but also hope to hear from @brorson on this and on calculation of the number of fourier terms to take as a function of m and q.

[brorson]

but also hope to hear from @brorson on this and on calculation of the number of fourier terms to take as a function of m and q.

Done. If I can do anything to help with the review process or merge let me know. Once you have a stable branch please tell me and I can check it out and run some of my other tests on it.

[steppi]

but also hope to hear from @brorson on this and on calculation of the number of fourier terms to take as a function of m and q.

Done. If I can do anything to help with the review process or merge let me know. Once you have a stable branch please tell me and I can check it out and run some of my other tests on it.

Thanks @brorson! I think this is getting into good shape now. I'll ping you again if we need more help. I will also run out some more tests against Mathematica's implementation.

[steppi]

@fbourgey, I think that should address everything from your first pass.

[fbourgey]

Minor things

[brorson]

Good by me.

Stuart

[steppi]

@fbourgey, can you give this one last look?

[fbourgey]

This LGTM, modulo with last comments.

[fbourgey]

Do we need a test for AngleUnitPolicy::Degrees as well?

[steppi]

Yeah, we can add some if you like.

[fbourgey]

what is the plan for this file? what about backward-compatibility?

should we not keep old tests for smoke coverage?

[steppi]

This file remains for now because it's still being used for the modified (radial) mathieu functions, and as a large order fallback for mathieu_cv. Once the modified mathieu functions are refactored like this, and we figure out something else for the fallback, then we can get rid of the file. We can't really keep the old tests because they rely on computing the functions from end to end and xsf doesn't offer that capability. The full implementations for CPU live in SciPy and rely on kernels from here that the CuPy implementations will also use. All of those old test cases (using degrees as angular unit) remain in SciPy.