Skip to content

Replace GSL quadrature and Levin sum in hyper1F1 with fortnum (2/5) - #30

Draft
krystophny wants to merge 5 commits into
mainfrom
drop-gsl-fortnum-hyper1f1
Draft

Replace GSL quadrature and Levin sum in hyper1F1 with fortnum (2/5)#30
krystophny wants to merge 5 commits into
mainfrom
drop-gsl-fortnum-hyper1f1

Conversation

@krystophny

@krystophny krystophny commented Jun 14, 2026

Copy link
Copy Markdown
Member

Merge order

Stack sequence: #29 -> #30 -> #31 -> #32 -> #33.

Each PR is now based on main individually, so its diff is cumulative against main and shares code with its stack neighbors. Merge in the sequence above. Once a predecessor merges into main, the next PR rebases onto it and its diff shrinks to its own increment.


Second in the fortnum migration stack (base: the Bessel PR). Moves the C confluent-hypergeometric evaluation off GSL.

Dependency removed: GSL adaptive quadrature and Levin u-sum acceleration in src/hyper1F1.c (gsl_integration_qag, gsl_integration_workspace, gsl_set_error_handler_off, gsl_sum_levin_u_*). The boundary integral now uses fortnum_integrate_qag with the Gauss-Kronrod 21 rule, and the Kummer-series acceleration uses fortnum_levin_u_accel. Both fortnum routines take the integrand plus an opaque context, so the existing quad_params struct is forwarded unchanged and no workspace is allocated. float.h is added for DBL_EPSILON, which was previously pulled in transitively through a GSL header.

Verification

Build the shared library at this tip:

$ cmake -S . -B build -G Ninja && cmake --build build --target mephit
[8/9] Building C object CMakeFiles/mephit.dir/src/hyper1F1.c.o
[9/9] Linking CXX shared library lib/libmephit.so

No GSL integration or sum symbols remain undefined; fortnum's routines are defined in the library:

$ nm -DuC build/lib/libmephit.so | grep -i 'gsl_integration_qag\|gsl_sum_levin\|gsl_integration_workspace'
(no output)
$ nm -DC build/lib/libmephit.so | grep -i 'fortnum_integrate_qag\|fortnum_levin_u'
0000000000115760 T fortnum_integrate_qag
0000000000115830 T fortnum_levin_u_accel

The only GSL symbols still referenced belong to the C call sites moved in PR 3 of this stack (glfixed, gsl_strerror). As in PR 1, the .x executables link only after PR 4 removes the empty FFTW link.

Note: fortnum pin updated to current main (92de6e9) after a fortnum history rewrite; old shas no longer resolve.

Drop the gsl_sf_bessel_In_array C binding in mephit_pert and compute the
modified Bessel functions through fortnum's bessel_in_array. fortnum fills
I_0..I_nmax in one pass, so the orders |m|-1, |m|, |m|+1 are recovered with
the symmetry I_{-n}(x) = I_n(x); for pol_mode = 0 the order -1 maps to I_1.

Wire fortnum via CMake FetchContent at the pinned tag, guarded so a build that
also pulls fortnum transitively through libneo does not declare the target
twice. Add the fortnum include directory for the C/C++ sources and link the
fortnum target. GSL stays linked until the remaining call sites move over.
@krystophny
krystophny marked this pull request as draft June 14, 2026 07:56
@krystophny
krystophny force-pushed the drop-gsl-fortnum-hyper1f1 branch from 8024bf9 to 12048cb Compare June 14, 2026 09:33
Route the a = 1 confluent hypergeometric evaluation through the fortnum C ABI:
the boundary integral uses fortnum_integrate_qag with the Gauss-Kronrod 21
rule, and the Kummer-series acceleration uses fortnum_levin_u_accel. Both
fortnum routines take the integrand plus an opaque context, so the existing
quad_params struct is forwarded unchanged and no workspace is allocated. Add
float.h for DBL_EPSILON, which was previously pulled in transitively through a
GSL header.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant