zeros of the kelvin functions

build.rs

@@ -379,6 +379,13 @@ void fcszo(int kf, int nt, cdouble *zo) {
     xsf::fcszo(kf, nt, zo);
 }"#;
 
+// kelvin.h
+
+const _CPP_KLVNZO: &str = r#"
+void klvnzo(int nt, int kd, double *zo) {
+    xsf::klvnzo(nt, kd, zo);
+}"#;
+
 // legendre.h
 
 const _CPP_ASSOC_LEGENDRE_P: &str = r#"
@@ -462,6 +469,10 @@ const WRAPPER_SPECS_CUSTOM: &[WrapperSpecCustom] = &[
         pattern: r"fcszo",
         cpp: _CPP_FCSZO,
     },
+    WrapperSpecCustom {
+        pattern: r"klvnzo",
+        cpp: _CPP_KLVNZO,
+    },
     WrapperSpecCustom {
         pattern: r"assoc_legendre_p_(0|1)",
         cpp: _CPP_ASSOC_LEGENDRE_P,

src/kelvin.rs

@@ -1,5 +1,7 @@
 use crate::bindings;
 use crate::bindings::xsf_impl;
+use alloc::vec::Vec;
+use core::ffi::c_int;
 use num_complex::Complex;
 
 xsf_impl!(ber, (x: f64), "Kelvin function `ber`");
@@ -14,23 +16,106 @@ xsf_impl!(keip, (x: f64), "Derivative of the Kelvin function `kei`");
 
 /// Kelvin functions as complex numbers
 ///
-/// # Arguments
-/// - `x` - Real argument
-///
-/// # Returns
-/// - *Be*: Value of the Kelvin function [`ber`] + *i* [`bei`]
-/// - *Ke*: Value of the Kelvin function [`ker`] + *i* [`kei`]
-/// - *Be*': Derivative of the Kelvin function [`berp`] + *i* [`beip`]
-/// - *Ke*': Derivative of the Kelvin function [`kerp`] + *i* [`keip`]
-pub fn kelvin(x: f64) -> (Complex<f64>, Complex<f64>, Complex<f64>, Complex<f64>) {
+/// Returned in a length-4 array of complex numbers:
+/// - *Be*: Value of [`ber`] + *i* [`bei`]
+/// - *Ke*: Value of [`ker`] + *i* [`kei`]
+/// - *Be*': Derivative of [`berp`] + *i* [`beip`]
+/// - *Ke*': Derivative of [`kerp`] + *i* [`keip`]
+pub fn kelvin(x: f64) -> [Complex<f64>; 4] {
     let mut be = bindings::complex_nan();
     let mut ke = bindings::complex_nan();
     let mut bep = bindings::complex_nan();
     let mut kep = bindings::complex_nan();
     unsafe {
         bindings::kelvin(x, &mut be, &mut ke, &mut bep, &mut kep);
     }
-    (be.into(), ke.into(), bep.into(), kep.into())
+    [be.into(), ke.into(), bep.into(), kep.into()]
+}
+
+enum KelvinFunction {
+    Ber = 1,
+    Bei = 2,
+    Ker = 3,
+    Kei = 4,
+    Berp = 5,
+    Beip = 6,
+    Kerp = 7,
+    Keip = 8,
+}
+
+#[inline(always)]
+fn klvnzo(nt: usize, kd: KelvinFunction) -> Vec<f64> {
+    assert!(nt <= c_int::MAX as usize);
+
+    let mut zs = alloc::vec![0.0; nt];
+    unsafe {
+        bindings::klvnzo(nt as c_int, kd as c_int, zs.as_mut_ptr());
+    }
+    zs
+}
+
+/// First `nt` zeros of Kelvin function [`ber`]
+pub fn ber_zeros(nt: usize) -> Vec<f64> {
+    klvnzo(nt, KelvinFunction::Ber)
+}
+
+/// First `nt` zeros of Kelvin function [`bei`]
+pub fn bei_zeros(nt: usize) -> Vec<f64> {
+    klvnzo(nt, KelvinFunction::Bei)
+}
+
+/// First `nt` zeros of Kelvin function [`ker`]
+pub fn ker_zeros(nt: usize) -> Vec<f64> {
+    klvnzo(nt, KelvinFunction::Ker)
+}
+
+/// First `nt` zeros of Kelvin function [`kei`]
+pub fn kei_zeros(nt: usize) -> Vec<f64> {
+    klvnzo(nt, KelvinFunction::Kei)
+}
+
+/// First `nt` zeros of Kelvin function derivative [`berp`]
+pub fn berp_zeros(nt: usize) -> Vec<f64> {
+    klvnzo(nt, KelvinFunction::Berp)
+}
+
+/// First `nt` zeros of Kelvin function derivative [`beip`]
+pub fn beip_zeros(nt: usize) -> Vec<f64> {
+    klvnzo(nt, KelvinFunction::Beip)
+}
+
+/// First `nt` zeros of Kelvin function derivative [`kerp`]
+pub fn kerp_zeros(nt: usize) -> Vec<f64> {
+    klvnzo(nt, KelvinFunction::Kerp)
+}
+
+/// First `nt` zeros of Kelvin function derivative [`keip`]
+pub fn keip_zeros(nt: usize) -> Vec<f64> {
+    klvnzo(nt, KelvinFunction::Keip)
+}
+
+/// First `nt` zeros of all Kelvin functions and their derivatives
+///
+/// Returned in a length-8 array of vectors of length nt. The array contains the vectors of zeros of
+/// - [`ber`]
+/// - [`bei`]
+/// - [`ker`]
+/// - [`kei`]
+/// - [`berp`]
+/// - [`beip`]
+/// - [`kerp`]
+/// - [`keip`]
+pub fn kelvin_zeros(nt: usize) -> [Vec<f64>; 8] {
+    [
+        klvnzo(nt, KelvinFunction::Ber),
+        klvnzo(nt, KelvinFunction::Bei),
+        klvnzo(nt, KelvinFunction::Ker),
+        klvnzo(nt, KelvinFunction::Kei),
+        klvnzo(nt, KelvinFunction::Berp),
+        klvnzo(nt, KelvinFunction::Beip),
+        klvnzo(nt, KelvinFunction::Kerp),
+        klvnzo(nt, KelvinFunction::Keip),
+    ]
 }
 
 #[cfg(test)]
@@ -84,7 +169,167 @@ mod tests {
         testing::test::<(Complex<f64>, Complex<f64>, Complex<f64>, Complex<f64>), _>(
             "kelvin",
             "d-cd_cd_cd_cd",
-            |x: &[f64]| kelvin(x[0]),
+            |x: &[f64]| kelvin(x[0]).into(),
+        );
+    }
+
+    // Kelvin function zeros tests
+
+    fn assert_allclose(actual: &[f64], expected: &[f64], atol: f64) {
+        assert_eq!(actual.len(), expected.len());
+        for (&a, &e) in actual.iter().zip(expected.iter()) {
+            assert!(
+                (a - e).abs() <= atol,
+                "actual: {}, expected: {}, diff: {}",
+                a,
+                e,
+                (a - e).abs()
+            );
+        }
+    }
+
+    /// Ported from `scipy.special.tests.test_basic.TestKelvin.test_ber_zeros`
+    #[test]
+    fn test_ber_zeros() {
+        let ber = ber_zeros(5);
+        assert_allclose(
+            &ber,
+            &[2.84892, 7.23883, 11.67396, 16.11356, 20.55463],
+            1.5e-4,
+        );
+    }
+
+    /// Ported from `scipy.special.tests.test_basic.TestKelvin.test_bei_zeros`
+    #[test]
+    fn test_bei_zeros() {
+        // Abramowitz & Stegun, Table 9.12
+        let bi = bei_zeros(5);
+        assert_allclose(
+            &bi,
+            &[5.02622, 9.45541, 13.89349, 18.33398, 22.77544],
+            1.5e-4,
+        );
+    }
+
+    /// Ported from `scipy.special.tests.test_basic.TestKelvin.test_ker_zeros`
+    #[test]
+    fn test_ker_zeros() {
+        let ker = ker_zeros(5);
+        assert_allclose(
+            &ker,
+            &[1.71854, 6.12728, 10.56294, 15.00269, 19.44381],
+            1.5e-4,
+        );
+    }
+
+    /// Ported from `scipy.special.tests.test_basic.TestKelvin.test_kei_zeros`
+    #[test]
+    fn test_kei_zeros() {
+        let kei = kei_zeros(5);
+        assert_allclose(
+            &kei,
+            &[3.91467, 8.34422, 12.78256, 17.22314, 21.66464],
+            1.5e-4,
+        );
+    }
+
+    /// Ported from `scipy.special.tests.test_basic.TestKelvin.test_berp_zeros`
+    #[test]
+    fn test_berp_zeros() {
+        let brp = berp_zeros(5);
+        assert_allclose(
+            &brp,
+            &[6.03871, 10.51364, 14.96844, 19.41758, 23.86430],
+            1.5e-4,
+        );
+    }
+
+    /// Ported from `scipy.special.tests.test_basic.TestKelvin.test_beip_zeros`
+    #[test]
+    fn test_beip_zeros() {
+        let bip = beip_zeros(5);
+        assert_allclose(
+            &bip,
+            &[
+                3.772_673_304_934_953,
+                8.280_987_849_760_042,
+                12.742_147_523_633_703,
+                17.193_431_752_512_54,
+                21.641_143_941_167_325,
+            ],
+            1.5e-8,
+        );
+    }
+
+    /// Ported from `scipy.special.tests.test_basic.TestKelvin.test_kerp_zeros`
+    #[test]
+    fn test_kerp_zeros() {
+        let kerp = kerp_zeros(5);
+        assert_allclose(
+            &kerp,
+            &[2.66584, 7.17212, 11.63218, 16.08312, 20.53068],
+            1.5e-4,
+        );
+    }
+
+    /// Ported from `scipy.special.tests.test_basic.TestKelvin.test_keip_zeros`
+    #[test]
+    fn test_keip_zeros() {
+        let keip = keip_zeros(5);
+        assert_allclose(
+            &keip,
+            &[4.93181, 9.40405, 13.85827, 18.30717, 22.75379],
+            1.5e-4,
+        );
+    }
+
+    /// Ported from `scipy.special.tests.test_basic.TestKelvin.test_kelvin_zeros`
+    #[test]
+    fn test_kelvin_zeros() {
+        // numbers come from 9.9 of A&S pg. 381
+        let tmp = kelvin_zeros(5);
+        let [berz, beiz, kerz, keiz, berpz, beipz, kerpz, keipz] = tmp;
+
+        assert_allclose(
+            &berz,
+            &[2.84892, 7.23883, 11.67396, 16.11356, 20.55463],
+            1.5e-4,
+        );
+        assert_allclose(
+            &beiz,
+            &[5.02622, 9.45541, 13.89349, 18.33398, 22.77544],
+            1.5e-4,
+        );
+        assert_allclose(
+            &kerz,
+            &[1.71854, 6.12728, 10.56294, 15.00269, 19.44382],
+            1.5e-4,
+        );
+        assert_allclose(
+            &keiz,
+            &[3.91467, 8.34422, 12.78256, 17.22314, 21.66464],
+            1.5e-4,
+        );
+        assert_allclose(
+            &berpz,
+            &[6.03871, 10.51364, 14.96844, 19.41758, 23.86430],
+            1.5e-4,
+        );
+        assert_allclose(
+            &beipz,
+            // table from 1927 had 3.77320 but this is more accurate
+            &[3.77267, 8.28099, 12.74215, 17.19343, 21.64114],
+            1.5e-4,
+        );
+        assert_allclose(
+            &kerpz,
+            &[2.66584, 7.17212, 11.63218, 16.08312, 20.53068],
+            1.5e-4,
+        );
+        assert_allclose(
+            &keipz,
+            &[4.93181, 9.40405, 13.85827, 18.30717, 22.75379],
+            1.5e-4,
         );
     }
 }

src/lib.rs

@@ -83,7 +83,10 @@ pub use iv_ratio::{iv_ratio, iv_ratio_c};
 
 // kelvin.h
 mod kelvin;
-pub use kelvin::{bei, beip, ber, berp, kei, keip, kelvin, ker, kerp};
+pub use kelvin::{
+    bei, bei_zeros, beip, beip_zeros, ber, ber_zeros, berp, berp_zeros, kei, kei_zeros, keip,
+    keip_zeros, kelvin, kelvin_zeros, ker, ker_zeros, kerp, kerp_zeros,
+};
 
 // lambertw.h
 mod lambertw;