Review request -- Mathieu fcns by brorson · Pull Request #60 · scipy/xsf

brorson · 2025-09-06T01:59:01Z

This is a request to try the new Mathieu fcn impl. This new fcn impl has been discussed in #47 . I have updated the docstrings and other comment in this checkin, but the functionality has otherwise remained the same as that shown in the bug report.

I have run these fcns against the Scipy test suite using "spin test" and the test suite runs to completion. The Mathieu fcns fail, but that is expected since this is a new impl and the input arg is in radians, not degrees.

Please review and comment. My goal here is to request a review and find out what I need to change.

Stuart

lucascolley · 2025-09-06T11:44:33Z

I'm afraid I can't help with C++ expertise, but I wonder whether you could open a draft PR to SciPy which shows the tests running (albeit with failures)? No rush, just thought it might be worth it whilst you are waiting for review here and have some momentum. Let me know if you would like help with making the necessary changes to the SciPy submodule etc.

brorson · 2025-09-07T23:37:52Z

Thank you for your suggestion, @lucascolley. I will do that.

Meanwhile, I see the CI thing is running on my PR and I have build failures. I fixed one -- I was using the suffix 'd' when declaring double literals which turns out to be non-standard. I fixed that one and pushed to my forked repo. Not sure if CI runs continuously but I will continue fixing any bugs it finds.

Stuart

lucascolley · 2025-09-08T06:48:04Z

It is a GitHub limitation that since this is your first PR, CI runs require approval from someone with write access on the repo.

brorson · 2025-09-10T12:45:00Z

OK, I just submitted a PR on the main Scipy branch. The pull request includes:

My changes to the meson.build file which are required to successfully build the Mathieu stuff. The change is imply to add the Lapack dependencies to ufuncs and other build targets per @steppi 's suggestion.
I also mentioned the failing tests. The Mathieu tests now fail since the implementation is different and returns different results from the old impl. (The old tests didn't check correctness; they just checked that the output of the fcns didn't change.) The required action item here is to bypass the tests for the time being so the impls can be integrated into the codebase. I also offer to provide some Python-based tests written by my students which can be used in the future. Maybe it would make sense to have a Scipy branch we can work in and not bother the trunk?

I also submitted a new PR on the xsf repo to pick up a bugfix I implemented after looking at some of my last CI test failures: I had suffixed some double literals using 'd', like this: 2.0d. This is apparently not standard C++ so made some of the compilers barf. My mistake -- I fixed it and issued a new PR to you.

Sorry for the multiple PRs -- I'm new to the Scipy dev chain and so will fumble around a little bit as I try to get this change into the codebase.

Stuart

steppi · 2025-09-10T15:08:47Z

Thanks @brorson. I plan to spend some time in the next few weeks to improve testing, and we can add some test cases based on references from https://github.com/mathieuandspheroidalwavefunctions computed with quad precision. In the meantime, your main.c should be moved into the tests folder and turned into a proper Catch2 test file named something like xsf/tests/miscellaneous_tests/test_mathieu_funcs.cpp that will get picked up by the test harness. You don't have to use parquet tables, you can just hardcode the golden values into the test files. Let me know if you need any help with this. If you do this, I can set the tolerances in the relevant tests in xsf/tests/scipy_special_tests to infinity until reference values are recomputed.

rgommers · 2025-09-10T15:41:59Z

It is a GitHub limitation that since this is your first PR, CI runs require approval from someone with write access on the repo.

The usual workaround is to open a separate PR fixing a typo or other such trivial thing somewhere in the repo; once that's merged, CI will run here.

brorson · 2025-09-10T22:47:53Z

Thank you, @steppi. Regarding the tests in main.c, sure I can move them to a file called "test_mathieu_fcns.cpp" and move them to xsf/tests/miscellaneous_tests/. Give me a day or two.

Once I do that I will check the mods into my GitHub repo and issue another PR so you know they're there.

brorson · 2025-09-10T23:57:26Z

@steppi, per your comments I moved the tests from main.c to a file xsf/tests/scipy_special_tests/test_mathieu.cpp in my local directory. But now I am looking at the other cpp files in that dir. It's clear that I need to do some homework and understand the formatting of the test files and how they are processed and then executed. Meanwhile, can I ask you to point me to some test files which are close to my use-case? My near-term goal is to just have a file which is picked up by the CI thing, runs some functions which test against golden values hard-coded into the file, then sets the right variables to tell the CI thing about passes/failures. Is there a particular file or files I can look at as an example?

Thanks,
Stuart

steppi · 2025-09-11T00:51:44Z

@steppi, per your comments I moved the tests from main.c to a file xsf/tests/scipy_special_tests/test_mathieu.cpp in my local directory. But now I am looking at the other cpp files in that dir. It's clear that I need to do some homework and understand the formatting of the test files and how they are processed and then executed. Meanwhile, can I ask you to point me to some test files which are close to my use-case? My near-term goal is to just have a file which is picked up by the CI thing, runs some functions which test against golden values hard-coded into the file, then sets the right variables to tell the CI thing about passes/failures. Is there a particular file or files I can look at as an example?

Thanks, Stuart

Thanks. I should have some time tomorrow to write a quick template for how the test file could be written and post it here

main.c to test_mathieu.cpp in the tests dir. I ran clang-format on all .h files.

brorson · 2025-09-11T01:07:52Z

Thank you.

jorenham

This looks pretty impressive!

I'm by no means an expert on the underlying maths or C++, but I left some (superficial) comments anyway; I hope you don't mind. I wasn't able to look at everything though, given the size and all.

One thing general pattern I've seen, is that there's quite a bit of code duplication. For reviewers that's not a lot of fun. So if you want to help speed up the review process, we'd appreciate it if you could try do something about the duplication. And there's no need to go wild with macro's; a couple of helper functions here and there could already help a lot.

I've also seen a several malloc and free calls. But since we're in C++ land, I don't think we need those. There's even the fancy std::mdspan you could use for vector- and matrix stuff, since we bundle a backport of it in kokkos/mdspan.hpp, which I imagine could be very useful here.

jorenham · 2025-09-19T20:45:11Z

+template <typename T>
+T cem_cva(T m, T q) {


since it internally casts to- and from double, I don't think that there's any advantage to having it generic like this

You are right, but this is how it was done in the old code, so I am reluctant to change it. This fcn is just a wrapper -- an adaption layer to the rest of Scipy.

jorenham · 2025-09-19T20:57:07Z

+    /*-----------------------------------------------
+    This creates the recurrence relation matrix for
+    the even-even Mathieu fcns (ce_2n).
+    Inputs:
+    N = matrix size (related to max order desired).
+    q = shape parameter.
+    Output:
+    A = recurrence matrix (must be calloc'ed in caller).
+    Return:
+    return code = 0 if OK.
+    -------------------------------------------------*/


nitpick: I think it would a bit easier to read if it would be a bit wider, so that there aren't as many linebreaks needed.

True, but my screen-size is limited since I work on a laptop and not a nice 4-screen work area. I prefer to keep the comment as is, if that's OK.

jorenham · 2025-09-19T21:59:28Z

+            *a = std::numeric_limits<double>::quiet_NaN();
+            return SF_ERROR_DOMAIN;


it's usually better to fail loudly in case of invalid input

This tries to replicate how the existing impl works. I am trying to replace the existing impl without too many surprises.

jorenham · 2025-09-19T22:04:05Z

+namespace mathieu {
+
+    //==================================================================
+    int mathieu_ce(int m, double q, double v, double *ce, double *ced) {


xsf uses pass-by-reference in case of multiple scalar outputs

Yes. However, this is in my internal code, and the old code works the same way. The adaption layer in xsf/mathieu.h adapts my internal stuff to the rest of xsf.

brorson · 2025-09-20T11:40:54Z

Thank you for the suggestions, @jorenham . I will go down the list and fix up my code based on your comments.

Regarding code duplication, yes I thought of using macros for the repeated stuff but since the repeated stuff is multi-line I figured any macro would be multi-line and hard to manage. However, I'll think about how to use macros again based on your comments, and may also see if I can stick some of the stuff into a separate inline fcn. I will want to comment heavily any macro or fcn to make clear what the calc is trying to do.

Regarding malloc() & free(), yeah, I am more of a C programmer than a C++ programmer. Does Scipy allow use of smart pointers? I see auto_ptr has been removed from C++ ... is there a preferred alternative?

(I will browse the codebase to see what is used, but it's always easier to ask ...)

Regarding failing loudly, do you mean throw an exception? I may be wrong, but since the xsf stuff supplies return codes I thought that was the preferred way to deal with errors. Please correct me if wrong.

Thx,
Stuart

brorson · 2025-09-20T11:54:51Z

Regarding code structure, mathieu.h living in the xsf directory provides the interface to the rest of Scipy. It is just a wrapper around my impls. My impls live in the mathieu subdirectory. Passing error codes within my stuff I have handled by just using integer flags, but I should fix that to conform to Scipy's error codes. Error codes passed back to Scipy from the interface fcns in mathieu.h are currently done using Scipy/xsf's predefined codes in error.h. I made a guess about what the different errors in the enum mean ... I hope I was close enough to be correct.

brorson · 2025-09-20T11:59:10Z

Regarding status of this project, I am currently bogged down with teaching and other responsibilities but hope to get back to this in a couple of weeks. Besides some clean-up per your suggestions (and my own list of things to clean up) I need to follow @rgommers suggestion about how to create a pull request which makes changes to two repos -- the main Scipy repo (where the Meson build config lives) and the xsf repo. I need to declare some dependencies in the build config stuff. That's waiting for me to do.

jorenham · 2025-09-23T11:11:10Z

Does Scipy allow use of smart pointers? I see auto_ptr has been removed from C++ ... is there a preferred alternative?

The cases I've seen could be replaced with a std::vector (docs) that you can pass around using .data() (docs)

I may be wrong, but since the xsf stuff supplies return codes I thought that was the preferred way to deal with errors

Anything user-facing that receives incorrect input returns NaN, generally speaking.
My reasoning here is that the major downside of return codes is that they're very easy to (accidentally) ignore, which is difficult for debugging. Exceptions must be explicitly dealt with, and if not, then you'll notice immediately, making it far less likely to miss things.
So since we're mostly talking about non-user-facing stuff here (unless I'm missing something), I'd choose exceptions over return codes, if I were in your shoes.
But either way, I meant this as a loose suggestion, not as a hard requirement or something.

functions to replace duplicated code.

and that's what I have.

fix the failed builds in CI.

…e it.

…hing.

and lapack which are requred for Mathieu fcns. Also fix the same build problem when using pixi. Also rerun clang-format on all files.

as well as the entirety of scipy.

brorson · 2025-12-26T16:04:31Z

Just wanted to bump this up. I have made a number of changes to my Mathieu fcn implementation. I have updated the code according to the reviewer's comments and made some of my own improvements. I also updated the build machinery so it finds and links to Lapack correctly. The code now builds successfully in the CI thing. Some XSF tests still fail, but that's expected since the old parquet files were generated using the old impl -- the old impl took inputs in degrees (!?!?!?) whereas my impl takes inputs in rad (more natural).

I would be happy to work with whoever is interested to update the xsf tests in a way which maintains mathematical correctness.

This impl also integrates into the main scipy trunk. The scipy tests also have some failures after modification for rad inputs, but the failures are for small testing tolerances. Again, I would be happy to work with the core devs to work through the scipy tests.

I am again requesting a review for possible merge into the main xsf trunk.

eigendecomposition on the full NxN recursion matrix, the new dstevd is for tridiagonal, symmetric matrices, so it runs significantly faster.

matrices (this didn't cause a bug because the problem was mostly in variable naming), updates to the choice of c, removal of the c=0 branch for the modified fcns, and others.

setting c. Minor comments added to make_matrix_tridiagonal.h.

col (m-2)/2 in the oe coeff fcn.

size scale with q somehow for better accuracy at high q. Also run clang-format on everything.

brorson · 2026-01-21T12:08:23Z

After some advice and some of my own work I think my reimplementation of the Mathieu functions can be pulled into xsf alone, not full Scipy.

The reason is that I made changes to the build system in xsf itself so the build stuff can find Lapack, which is required for the Mathieu build. Previously I needed to pull into Scipy to get Lapack.

I will close this pull request and open it under xsf instead.

brorson added 3 commits September 5, 2025 20:51

Update mathieu.h with improved docstrings and comments.

0adc5a6

Add new directory of Mathieu fcn impls.

9f62cf0

Make sure we return nan when error is encountered.

204e09d

Remove non-standard d in double constant literals.

b09c280

Merge branch 'scipy:main' into main

41c1f28

rgommers mentioned this pull request Sep 10, 2025

Add Lapack dependencies to meson.build for XSF scipy/scipy#23593

Closed

rgommers reviewed Sep 10, 2025

View reviewed changes

Comment thread include/xsf/mathieu/LICENSE Outdated

brorson added 2 commits September 10, 2025 19:28

Remove license file and place code under whatever license xsf uses.

93a41df

Merge branch 'main' of https://github.com/brorson/xsf_mathieu

b8911c1

I removed my Makefile since it's not used in the Scipy CI tool. I moved

5cbf45c

main.c to test_mathieu.cpp in the tests dir. I ran clang-format on all .h files.

Run clang-format on test_mathieu.cpp so it passes CI.

cdf598b

jorenham reviewed Sep 19, 2025

View reviewed changes

brorson added 4 commits September 24, 2025 18:59

Merge branch 'scipy:main' into main

cec5ca1

Replaced malloc/free with std::vector, updated some comments, added

9c850da

functions to replace duplicated code.

Merge branch 'main' of https://github.com/brorson/xsf_mathieu

72a5146

Merge branch 'scipy:main' into main

b25ba42

brorson added 15 commits November 13, 2025 14:26

Change dsyevd to dstevd since dstevd is for tridiagonal matrices,

5ddc975

and that's what I have.

Revert back to dsyevd. Use correct calling sig this time.

67314ef

Clean up files using clang-format before checking in.

c5e098b

Merge branch 'scipy:main' into main

7873a34

I used a more recent version of clang-format. Hopefully this will

78ce632

fix the failed builds in CI.

Modify meson.build so it includes Lapack as a dependency.

f3074d5

Don't use _Float128 if not defined on this platform.

536959f

Make _Float128 directive obey clang-format.

eac37fe

Fix some constant doubles -- don't use d suffix since Mac doesn't lik…

1d0ec45

…e it.

Run clang-format on mathieu_fcns to satisfy Scipy ci thing.

7eb1c86

Merge branch 'scipy:main' into main

d200bd4

Remove d suffix from constant doubles in test fcn to appease the CI t…

eb44cb3

…hing.

Fix meson build problem where standalone build of xsf doesn't find blas

8a401b3

and lapack which are requred for Mathieu fcns. Also fix the same build problem when using pixi. Also rerun clang-format on all files.

Updates to the build files again so I can build both stand-alone xsf

678e532

as well as the entirety of scipy.

Updates to mathieu.h and make_matrix.h based on reviewer comments.

2c0710f

brorson added 10 commits January 3, 2026 17:29

Change impl to use dstevd instead of dsyevd. The old dsyevd did an

ce4baf2

eigendecomposition on the full NxN recursion matrix, the new dstevd is for tridiagonal, symmetric matrices, so it runs significantly faster.

This creates the tridiagonal recursion matrix suitable for dstevd.

cb4b6d6

Udpate the tests to match the new API for the dstevd impl.

f08887b

I ran clang-format on the .h files since I forgot to do this last time.

caab938

Big changes including fixing row/col confusion in the eigenvector

3c17ce4

matrices (this didn't cause a bug because the problem was mostly in variable naming), updates to the choice of c, removal of the c=0 branch for the modified fcns, and others.

Run clang-format on *.h files.

36f71c9

Old s=0 code removed from mathieu_fcns. Also added s-max aproach to

80d930d

setting c. Minor comments added to make_matrix_tridiagonal.h.

Fix bug in se_2m+2 by changing coeff generation -- ask for matrix

d1da899

col (m-2)/2 in the oe coeff fcn.

Add get_partial_sum_N to mathieu_fcns. The goal is to make the matrix

d6b2c3f

size scale with q somehow for better accuracy at high q. Also run clang-format on everything.

Must update the mathieu test using clang-format.

bfd8ba5

brorson marked this pull request as draft January 21, 2026 11:46

brorson closed this Jan 21, 2026

rgommers mentioned this pull request Apr 26, 2026

Updated Mathieu function implementation #99

Closed

		*a = std::numeric_limits<double>::quiet_NaN();
		return SF_ERROR_DOMAIN;

Uh oh!

Uh oh!

Conversation

brorson commented Sep 6, 2025

Uh oh!

lucascolley commented Sep 6, 2025

Uh oh!

brorson commented Sep 7, 2025

Uh oh!

lucascolley commented Sep 8, 2025

Uh oh!

brorson commented Sep 10, 2025

Uh oh!

steppi commented Sep 10, 2025

Uh oh!

rgommers commented Sep 10, 2025

Uh oh!

Uh oh!

brorson commented Sep 10, 2025

Uh oh!

brorson commented Sep 10, 2025

Uh oh!

steppi commented Sep 11, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

brorson commented Sep 11, 2025

Uh oh!

jorenham left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

brorson commented Sep 20, 2025

Uh oh!

brorson commented Sep 20, 2025

Uh oh!

brorson commented Sep 20, 2025

Uh oh!

jorenham commented Sep 23, 2025

Uh oh!

brorson commented Dec 26, 2025

Uh oh!

brorson commented Jan 21, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

steppi commented Sep 11, 2025 •

edited

Loading