diff --git a/build.rs b/build.rs index 170e3ae..3544f1e 100644 --- a/build.rs +++ b/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"), diff --git a/src/fresnel.rs b/src/fresnel.rs new file mode 100644 index 0000000..a6e4327 --- /dev/null +++ b/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 {} +} + +pub trait FresnelArg: sealed::Sealed { + type Output; + fn fresnel(self) -> (Self::Output, Self::Output); +} + +#[inline(always)] +fn c_c64_nan() -> bindings::root::std::complex { + 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 { + type Output = Complex; + + #[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(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, Complex) { + 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, Complex) { + 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()) +} diff --git a/src/lib.rs b/src/lib.rs index e4cbcda..9762a24 100644 --- a/src/lib.rs +++ b/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}; diff --git a/tests/test_functions.rs b/tests/test_functions.rs index 55ab717..17609e2 100644 --- a/tests/test_functions.rs +++ b/tests/test_functions.rs @@ -59,6 +59,28 @@ mod xsref { } } + impl TestOutput for (f64, f64) { + fn from_parquet_rows(rows: Vec>) -> Vec { + 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 { rows.into_iter() @@ -88,6 +110,33 @@ mod xsref { } } + impl TestOutput for (Complex, Complex) { + fn from_parquet_rows(rows: Vec>) -> Vec { + 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, Complex, Complex, Complex) { fn from_parquet_rows(rows: Vec>) -> Vec { 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, Complex)) => { + paste::paste! { + #[test] + fn $test_name() { + xsref::test::<(num_complex::Complex, num_complex::Complex), _>(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!( + [], + $f, + "d-d_d", + |x: &[f64]| xsf::$f(x[0]), + (f64, f64) + ); + } + }; + (@single $f:ident, "D->DD") => { + paste::paste! { + _test!( + [], + $f, + "cd-cd_cd", + |x: &[f64]| xsf::$f(num_complex::Complex::new(x[0], x[1])), + (Complex, Complex) + ); + } + }; + (@single $f:ident, "d->DD") => { + paste::paste! { + _test!( + [], + $f, + "d-cd_cd", + |x: &[f64]| xsf::$f(x[0]), + (Complex, Complex) + ); + } + }; (@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");