Skip to content

Conversation

@WarrenWeckesser
Copy link
Member

Closes gh-46.

Thanks to github user hchau630 for identifying the problem and suggesting the fix.

Closes scipygh-46.

Thanks to github user hchau630 for identifying the problem
and suggesting the fix.
@WarrenWeckesser
Copy link
Member Author

WarrenWeckesser commented Oct 29, 2025

I didn't add new unit tests. At the moment, I'm not sure what the recommended way to do that is here in xsf. Edit: tests added.

In the mean time, here's a C++ program that accepts the input parameters of cyl_bessel_k() on the command line and prints the result. It allows a quick comparison to the result computed with mpmath.

cyl_bessel_k_demo.cpp

// Requires at least C++17.

#include <complex>
#include <iomanip>
#include <iostream>
#include <xsf/bessel.h>

int main(int argc, char *argv[])
{
    using namespace std;

    double nu = 0.0;
    double re = 0.0;
    double im = 0.0;

    if (argc > 1) {
        nu = stod(argv[1]);
        if (argc > 2) {
            re = stod(argv[2]);
            if (argc > 3) {
                im = stod(argv[3]);
            }
        }
    }
    complex<double> z(re, im);
    complex<double> w = xsf::cyl_bessel_k(nu, z);

    cout << setprecision(15);
    cout << w.real() << " " << w.imag() << endl;
    return 0;
}

For example, with the changes in this PR:

% ./cyl_bessel_k_demo 0 680 680 
-4.55373003294427e-298 -1.87872701010989e-297
% python3 -c "from mpmath import mp; mp.dps = 100; z = 680+680j; w = complex(mp.besselk(0, z)); print(f'{w.real:.15e} {w.imag:.15e}')"
-4.553730032944803e-298 -1.878727010109855e-297
% 
% ./cyl_bessel_k_demo 0 680 1000
1.90168487199878e-298 -1.71341234147958e-297
% python3 -c "from mpmath import mp; mp.dps = 100; z = 680+1000j; w = complex(mp.besselk(0, z)); print(f'{w.real:.15e} {w.imag:.15e}')"
1.901684871999608e-298 -1.713412341479591e-297

@WarrenWeckesser
Copy link
Member Author

After poking around some more, I see there are tests that don't rely on those parquet files. For example, @jorenham's pull request added test_lqmn.cpp. So I can put together a simple test for this change, and we can iterate on the preferred test idioms to use.

@WarrenWeckesser
Copy link
Member Author

I added unit tests. Since these are new tests that are not related to old scipy tests, I put the test file in a new subdirectory, xsf_tests (instead of scipy_special_tests).

@WarrenWeckesser
Copy link
Member Author

Windows failures are not new: #70.
The GPU job will be canceled in about 24 hours. ☹️

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

bug Something isn't working

Projects

None yet

Development

Successfully merging this pull request may close these issues.

BUG: Incorrect translation of AMOS Fortran code resulting in incorrect/NaN outputs

2 participants