erf_zeros

build.rs

@@ -18,6 +18,7 @@ const WRAPPER_INCLUDES: &[&str] = &[
     "cephes/round.h",
     "cephes/spence.h",
     "cephes/unity.h",
+    "specfun/specfun.h",
     "airy.h",
     "alg.h",
     "bessel.h",
@@ -430,6 +431,17 @@ double lgam1p(double x) {
     return xsf::cephes::lgam1p(x);
 }"#;
 
+// specfun/specfun.h
+
+const _CPP_SPECFUN_CERZO: &str = r#"
+void cerzo(int nt, cdouble *zo) {
+    std::vector<std::complex<double>> czo(nt);
+    xsf::specfun::cerzo(nt, czo.data());
+    for (int i = 0; i < nt; i++) {
+        zo[i] = cdouble(czo[i]);
+    }
+}"#;
+
 // airy.h
 
 const _CPP_AIRYZO: &str = r#"
@@ -628,6 +640,10 @@ const WRAPPER_SPECS_CUSTOM: &[WrapperSpecCustom] = &[
         pattern: r"lgam1p",
         cpp: _CPP_CEPHES_UNITY,
     },
+    WrapperSpecCustom {
+        pattern: r"cerzo",
+        cpp: _CPP_SPECFUN_CERZO,
+    },
     WrapperSpecCustom {
         pattern: r"airyzo",
         cpp: _CPP_AIRYZO,

src/lib.rs

@@ -248,6 +248,7 @@
 //! | [`erfi`]                   | Imaginary error function *-i erf(i z)*       |
 //! | [`erfinv`]                 | Inverse of [`erf`]                           |
 //! | [`erfcinv`]                | Inverse of [`erfc`]                          |
+//! | [`erf_zeros`]              | Zeros of [`erf`]                             |
 //! | [`wofz`]                   | Faddeeva function                            |
 //! | [`dawsn`]                  | Dawson's integral                            |
 //! | [`fresnel`]                | Fresnel integrals *S(z)* and *C(z)*          |

src/xsf/erf.rs

@@ -1,3 +1,4 @@
+use core::ffi::c_int;
 use num_complex::Complex;
 
 mod sealed {
@@ -80,6 +81,21 @@ pub fn erfi<T: ErfArg>(z: T) -> T {
     z.erfi()
 }
 
+/// Zeros of the error function `erf(z)` in the first quadrant of the complex plane
+///
+/// Zeros in the other quadrants can be obtained by using the symmetries *erf(-z) = erf(z)* and
+/// *erf(conj(z)) = conj(erf(z))*.
+///
+/// See [`erf`](crate::erf) for the error function itself.
+#[doc(alias = "cerzo")]
+pub fn erf_zeros<const NT: usize>() -> [Complex<f64>; NT] {
+    let mut zo = [f64::NAN.into(); NT];
+    unsafe {
+        crate::ffi::xsf::cerzo(NT as c_int, zo.as_mut_ptr());
+    }
+    zo.map(Into::into)
+}
+
 /// Dawson function `sqrt(pi)/2 * exp(-z^2) * erfi(z)` for real or complex input
 #[doc(alias = "dawson")]
 pub fn dawsn<T: ErfArg>(z: T) -> T {
@@ -99,83 +115,106 @@ pub fn voigt_profile(x: f64, sigma: f64, gamma: f64) -> f64 {
 
 #[cfg(test)]
 mod tests {
-    use super::*;
+    use crate::testing::np_assert_allclose;
     use crate::xsref;
-    use num_complex::{Complex, c64};
+    use num_complex::c64;
 
     // erf
 
     #[test]
     fn test_erf_f64() {
-        xsref::test::<f64, _>("erf", "d-d", |x: &[f64]| erf(x[0]));
+        xsref::test("erf", "d-d", |x: &[f64]| crate::erf(x[0]));
     }
 
     #[test]
     fn test_erf_c64() {
-        xsref::test::<Complex<f64>, _>("erf", "cd-cd", |x: &[f64]| erf(c64(x[0], x[1])));
+        xsref::test("erf", "cd-cd", |x: &[f64]| crate::erf(c64(x[0], x[1])));
     }
 
     // erfc
 
     #[test]
     fn test_erfc_f64() {
-        xsref::test::<f64, _>("erfc", "d-d", |x: &[f64]| erfc(x[0]));
+        xsref::test("erfc", "d-d", |x: &[f64]| crate::erfc(x[0]));
     }
 
     #[test]
     fn test_erfc_c64() {
-        xsref::test::<Complex<f64>, _>("erfc", "cd-cd", |x: &[f64]| erfc(c64(x[0], x[1])));
+        xsref::test("erfc", "cd-cd", |x: &[f64]| crate::erfc(c64(x[0], x[1])));
     }
 
     // erfcx
 
     #[test]
     fn test_erfcx_f64() {
-        xsref::test::<f64, _>("erfcx", "d-d", |x: &[f64]| erfcx(x[0]));
+        xsref::test("erfcx", "d-d", |x: &[f64]| crate::erfcx(x[0]));
     }
 
     #[test]
     fn test_erfcx_c64() {
-        xsref::test::<Complex<f64>, _>("erfcx", "cd-cd", |x: &[f64]| erfcx(c64(x[0], x[1])));
+        xsref::test("erfcx", "cd-cd", |x: &[f64]| crate::erfcx(c64(x[0], x[1])));
     }
 
     // erfi
 
     #[test]
     fn test_erfi_f64() {
-        xsref::test::<f64, _>("erfi", "d-d", |x: &[f64]| erfi(x[0]));
+        xsref::test("erfi", "d-d", |x: &[f64]| crate::erfi(x[0]));
     }
 
     #[test]
     fn test_erfi_c64() {
-        xsref::test::<Complex<f64>, _>("erfi", "cd-cd", |x: &[f64]| erfi(c64(x[0], x[1])));
+        xsref::test("erfi", "cd-cd", |x: &[f64]| crate::erfi(c64(x[0], x[1])));
+    }
+
+    // erf_zeros (from `scipy.speceial.tests.test_basic.TestErf.test_erf_zeros`)
+
+    #[test]
+    fn test_erf_zeros() {
+        // erz = special.erf_zeros(5)
+        let erz = crate::erf_zeros::<5>();
+        // erzr = array([1.45061616+1.88094300j,
+        //               2.24465928+2.61657514j,
+        //               2.83974105+3.17562810j,
+        //               3.33546074+3.64617438j,
+        //               3.76900557+4.06069723j])
+        let erzr = [
+            c64(1.45061616, 1.88094300),
+            c64(2.24465928, 2.61657514),
+            c64(2.83974105, 3.17562810),
+            c64(3.33546074, 3.64617438),
+            c64(3.76900557, 4.06069723),
+        ];
+        // assert_allclose(erz, erzr, atol=1.5e-4, rtol=0)
+        np_assert_allclose(&erz.map(|z| z.re), &erzr.map(|z| z.re), 0.0, 1.5e-4);
+        np_assert_allclose(&erz.map(|z| z.im), &erzr.map(|z| z.im), 0.0, 1.5e-4);
     }
 
     // dawsn
 
     #[test]
     fn test_dawsn_f64() {
-        xsref::test::<f64, _>("dawsn", "d-d", |x: &[f64]| dawsn(x[0]));
+        xsref::test("dawsn", "d-d", |x: &[f64]| crate::dawsn(x[0]));
     }
 
     #[test]
     fn test_dawsn_c64() {
-        xsref::test::<Complex<f64>, _>("dawsn", "cd-cd", |x: &[f64]| dawsn(c64(x[0], x[1])));
+        xsref::test("dawsn", "cd-cd", |x: &[f64]| crate::dawsn(c64(x[0], x[1])));
     }
 
     // wofz
 
     #[test]
     fn test_wofz() {
-        xsref::test::<Complex<f64>, _>("wofz", "cd-cd", |x: &[f64]| wofz(c64(x[0], x[1])));
+        xsref::test("wofz", "cd-cd", |x: &[f64]| crate::wofz(c64(x[0], x[1])));
     }
 
     // voigt_profile
 
     #[test]
     fn test_voigt_profile() {
-        xsref::test::<f64, _>("voigt_profile", "d_d_d-d", |x: &[f64]| {
-            voigt_profile(x[0], x[1], x[2])
+        xsref::test("voigt_profile", "d_d_d-d", |x: &[f64]| {
+            crate::voigt_profile(x[0], x[1], x[2])
         });
     }
 }

src/xsf/mod.rs

@@ -50,7 +50,7 @@ pub use binom::binom;
 pub use cdflib::gdtrib;
 pub use digamma::digamma;
 pub use ellip::{ellipe, ellipeinc, ellipj, ellipk, ellipkinc, ellipkm1};
-pub use erf::{dawsn, erf, erfc, erfcx, erfi, voigt_profile, wofz};
+pub use erf::{dawsn, erf, erf_zeros, erfc, erfcx, erfi, voigt_profile, wofz};
 pub use evalpoly::cevalpoly;
 pub use exp::{exp2, exp10, expm1};
 pub use expint::{exp1, expi, scaled_exp1};