fresnel integrals

build.rs

@@ -19,7 +19,7 @@ const XSF_HEADERS: &[&str] = &[
     "exp.h",
     "expint.h",
     // "fp_error_metrics.h",
-    // "fresnel.h",
+    "fresnel.h",
     "gamma.h",
     "hyp2f1.h",
     "iv_ratio.h",
@@ -133,7 +133,11 @@ const XSF_TYPES: &[(&str, &str)] = &[
     ("expi", "D->D"),
     ("scaled_exp1", "d->d"),
     // fresnel.h
-    //  TODO
+    //  TODO: `fcszo`: ii->[D]
+    ("fresnel", "d->dd"),
+    ("fresnel", "D->DD"),
+    ("modified_fresnel_plus", "d->DD"),
+    ("modified_fresnel_minus", "d->DD"),
     // gamma.h
     ("gamma", "d->d"),
     ("gamma", "D->D"),

src/fresnel.rs

@@ -0,0 +1,133 @@
+use crate::bindings;
+use num_complex::Complex;
+
+mod sealed {
+    use num_complex::Complex;
+
+    pub trait Sealed {}
+    impl Sealed for f64 {}
+    impl Sealed for Complex<f64> {}
+}
+
+pub trait FresnelArg: sealed::Sealed {
+    type Output;
+    fn fresnel(self) -> (Self::Output, Self::Output);
+}
+
+#[inline(always)]
+fn c_c64_nan() -> bindings::root::std::complex<f64> {
+    Complex::new(f64::NAN, f64::NAN).into()
+}
+
+impl FresnelArg for f64 {
+    type Output = f64;
+
+    #[inline(always)]
+    fn fresnel(self) -> (Self::Output, Self::Output) {
+        let mut fs = f64::NAN;
+        let mut fc = f64::NAN;
+        unsafe {
+            bindings::fresnel(self, &mut fs, &mut fc);
+        }
+        (fs, fc)
+    }
+}
+
+impl FresnelArg for Complex<f64> {
+    type Output = Complex<f64>;
+
+    #[inline(always)]
+    fn fresnel(self) -> (Self::Output, Self::Output) {
+        let mut fs = c_c64_nan();
+        let mut fc = c_c64_nan();
+
+        unsafe {
+            bindings::fresnel_1(self.into(), &mut fs, &mut fc);
+        }
+        (fs.into(), fc.into())
+    }
+}
+
+/// Fresnel integrals
+///
+/// The Fresnel integrals are defined as:
+///
+/// - S(z) = ∫₀ᶻ sin(π t² / 2) dt
+/// - C(z) = ∫₀ᶻ cos(π t² / 2) dt
+///
+/// where the integrals are taken from 0 to z.
+///
+/// # Arguments
+///
+/// - `z` - real- or complex-valued argument
+///
+/// # Returns
+///
+/// - `S` - S(z)
+/// - `C` - C(z)
+///
+/// # See also
+///
+/// - [`modified_fresnel_plus`] - Modified Fresnel positive integrals
+/// - [`modified_fresnel_minus`] - Modified Fresnel negative integrals
+pub fn fresnel<T: FresnelArg>(z: T) -> (T::Output, T::Output) {
+    z.fresnel()
+}
+
+/// Modified Fresnel positive integrals
+///
+/// Computes the modified Fresnel integrals:
+///
+/// - F₊(x) = ∫ₓ∞ exp(i t²) dt
+/// - K₊(x) = F₊(x) × exp(-i (x² + π/4)) / √π
+///
+/// # Arguments
+///
+/// - `x` - real-valued argument
+///
+/// # Returns
+///
+/// - `fp` - Integral F₊(x)
+/// - `kp` - Integral K₊(x)
+///
+/// # See also
+///
+/// - [`fresnel`] - Standard Fresnel integrals
+/// - [`modified_fresnel_minus`] - Modified Fresnel negative integrals
+pub fn modified_fresnel_plus(x: f64) -> (Complex<f64>, Complex<f64>) {
+    let mut fp = c_c64_nan();
+    let mut kp = c_c64_nan();
+    unsafe {
+        bindings::modified_fresnel_plus(x, &mut fp, &mut kp);
+    }
+    (fp.into(), kp.into())
+}
+
+/// Modified Fresnel negative integrals
+///
+/// Computes the modified Fresnel integrals:
+///
+/// - F₋(x) = ∫ₓ∞ exp(-i t²) dt
+/// - K₋(x) = F₋(x) × exp(i (x² + π/4)) / √π
+///
+/// # Arguments
+///
+/// - `x` - real-valued argument
+///
+/// # Returns
+///
+/// - `fm` - Integral F₋(x)
+/// - `km` - Integral K₋(x)
+///
+/// # See also
+///
+/// - [`fresnel`] - Standard Fresnel integrals
+/// - [`modified_fresnel_plus`] - Modified Fresnel positive integrals
+pub fn modified_fresnel_minus(x: f64) -> (Complex<f64>, Complex<f64>) {
+    let mut fm = c_c64_nan();
+    let mut km = c_c64_nan();
+    unsafe {
+        bindings::modified_fresnel_minus(x, &mut fm, &mut km);
+    }
+    (fm.into(), km.into())
+}

src/lib.rs

@@ -50,6 +50,10 @@ pub use exp::{exp2, exp10, expm1};
 mod expint;
 pub use expint::{exp1, expi, scaled_exp1};
 
+// fresnel.h
+mod fresnel;
+pub use fresnel::{fresnel, modified_fresnel_minus, modified_fresnel_plus};
+
 // gamma.h
 mod gamma;
 pub use gamma::{gamma, gammainc, gammaincc, gammainccinv, gammaincinv, gammaln, gammasgn};

tests/test_functions.rs

@@ -59,6 +59,28 @@ mod xsref {
         }
     }
 
+    impl TestOutput for (f64, f64) {
+        fn from_parquet_rows(rows: Vec<Vec<f64>>) -> Vec<Self> {
+            rows.into_iter().map(|row| (row[0], row[1])).collect()
+        }
+
+        fn error(actual: Self, expected: Self) -> f64 {
+            let errors = [
+                relative_error(actual.0, expected.0),
+                relative_error(actual.1, expected.1),
+            ];
+            errors.iter().fold(0.0, |acc, &x| acc.max(x))
+        }
+
+        fn magnitude(self) -> f64 {
+            (self.0.abs() + self.1.abs()) / 2.0
+        }
+
+        fn format_value(self) -> String {
+            format!("({:.6e}, {:.6e})", self.0, self.1)
+        }
+    }
+
     impl TestOutput for (f64, f64, f64, f64) {
         fn from_parquet_rows(rows: Vec<Vec<f64>>) -> Vec<Self> {
             rows.into_iter()
@@ -88,6 +110,33 @@ mod xsref {
         }
     }
 
+    impl TestOutput for (Complex<f64>, Complex<f64>) {
+        fn from_parquet_rows(rows: Vec<Vec<f64>>) -> Vec<Self> {
+            rows.into_iter()
+                .map(|row| (Complex::new(row[0], row[1]), Complex::new(row[2], row[3])))
+                .collect()
+        }
+
+        fn error(actual: Self, expected: Self) -> f64 {
+            let errors = [
+                complex_relative_error(actual.0, expected.0),
+                complex_relative_error(actual.1, expected.1),
+            ];
+            errors.iter().fold(0.0, |acc, &x| acc.max(x))
+        }
+
+        fn magnitude(self) -> f64 {
+            (self.0.norm() + self.1.norm()) / 2.0
+        }
+
+        fn format_value(self) -> String {
+            format!(
+                "({:.6e}+{:.6e}i, {:.6e}+{:.6e}i)",
+                self.0.re, self.0.im, self.1.re, self.1.im,
+            )
+        }
+    }
+
     impl TestOutput for (Complex<f64>, Complex<f64>, Complex<f64>, Complex<f64>) {
         fn from_parquet_rows(rows: Vec<Vec<f64>>) -> Vec<Self> {
             rows.into_iter()
@@ -422,6 +471,22 @@ macro_rules! _test {
             }
         }
     };
+    ($test_name:ident, $f:ident, $sig:literal, $test_fn:expr, (f64, f64)) => {
+        paste::paste! {
+            #[test]
+            fn $test_name() {
+                xsref::test::<(f64, f64), _>(stringify!($f), $sig, $test_fn);
+            }
+        }
+    };
+    ($test_name:ident, $f:ident, $sig:literal, $test_fn:expr, (Complex<f64>, Complex<f64>)) => {
+        paste::paste! {
+            #[test]
+            fn $test_name() {
+                xsref::test::<(num_complex::Complex<f64>, num_complex::Complex<f64>), _>(stringify!($f), $sig, $test_fn);
+            }
+        }
+    };
     ($test_name:ident, $f:ident, $sig:literal, $test_fn:expr, (f64, f64, f64, f64)) => {
         paste::paste! {
             #[test]
@@ -528,6 +593,39 @@ macro_rules! xsref_test {
             );
         }
     };
+    (@single $f:ident, "d->dd") => {
+        paste::paste! {
+            _test!(
+                [<test_ $f _ d>],
+                $f,
+                "d-d_d",
+                |x: &[f64]| xsf::$f(x[0]),
+                (f64, f64)
+            );
+        }
+    };
+    (@single $f:ident, "D->DD") => {
+        paste::paste! {
+            _test!(
+                [<test_ $f _ cd>],
+                $f,
+                "cd-cd_cd",
+                |x: &[f64]| xsf::$f(num_complex::Complex::new(x[0], x[1])),
+                (Complex<f64>, Complex<f64>)
+            );
+        }
+    };
+    (@single $f:ident, "d->DD") => {
+        paste::paste! {
+            _test!(
+                [<test_ $f _ cd>],
+                $f,
+                "d-cd_cd",
+                |x: &[f64]| xsf::$f(x[0]),
+                (Complex<f64>, Complex<f64>)
+            );
+        }
+    };
     (@single $f:ident, "d->dddd") => {
         paste::paste! {
             _test!(
@@ -713,6 +811,11 @@ xsref_test!(expi, "d->d", "D->D");
 xsref_test!(exp1, "d->d", "D->D");
 xsref_test!(scaled_exp1, "d->d");
 
+// fresnel.h
+xsref_test!(fresnel, "d->dd", "D->DD");
+xsref_test!(modified_fresnel_plus, "d->DD");
+xsref_test!(modified_fresnel_minus, "d->DD");
+
 // gamma.h
 xsref_test!(gamma, "d->d", "D->D");
 xsref_test!(gammainc, "dd->d");