diff --git a/src/cauchy_pv.rs b/src/cauchy_pv.rs new file mode 100644 index 0000000..c45ed71 --- /dev/null +++ b/src/cauchy_pv.rs @@ -0,0 +1,214 @@ +//! Cauchy principal value integration. +//! +//! Computes PV ∫ₐᵇ f(x)/(x-c) dx where c ∈ (a, b) using the subtraction +//! technique: rewrite as a regular integral plus an analytic logarithmic term. +//! +//! # Example +//! +//! ``` +//! use bilby::cauchy_pv::pv_integrate; +//! +//! // PV ∫₀¹ x²/(x - 0.3) dx = 0.8 + 0.09 * ln(7/3) +//! let exact = 0.8 + 0.09 * (7.0_f64 / 3.0).ln(); +//! let result = pv_integrate(|x| x * x, 0.0, 1.0, 0.3, 1e-10).unwrap(); +//! assert!((result.value - exact).abs() < 1e-8); +//! ``` + +use crate::adaptive::AdaptiveIntegrator; +use crate::error::QuadratureError; +use crate::result::QuadratureResult; + +/// Builder for Cauchy principal value integration. +/// +/// Computes PV ∫ₐᵇ f(x)/(x-c) dx via the subtraction technique: +/// +/// PV ∫ₐᵇ f(x)/(x-c) dx = ∫ₐᵇ \[f(x) - f(c)\]/(x-c) dx + f(c) · ln((b-c)/(c-a)) +/// +/// The first integral has a removable singularity at x = c and is computed +/// using adaptive quadrature with a break point at c. +#[derive(Debug, Clone)] +pub struct CauchyPV { + abs_tol: f64, + rel_tol: f64, + max_evals: usize, +} + +impl Default for CauchyPV { + fn default() -> Self { + Self { + abs_tol: 1.49e-8, + rel_tol: 1.49e-8, + max_evals: 500_000, + } + } +} + +impl CauchyPV { + /// Set absolute tolerance. + pub fn with_abs_tol(mut self, tol: f64) -> Self { + self.abs_tol = tol; + self + } + + /// Set relative tolerance. + pub fn with_rel_tol(mut self, tol: f64) -> Self { + self.rel_tol = tol; + self + } + + /// Set maximum number of function evaluations. + pub fn with_max_evals(mut self, n: usize) -> Self { + self.max_evals = n; + self + } + + /// Compute PV ∫ₐᵇ f(x)/(x-c) dx. + /// + /// The singularity c must be strictly inside (a, b). + pub fn integrate( + &self, + a: f64, + b: f64, + c: f64, + f: G, + ) -> Result, QuadratureError> + where + G: Fn(f64) -> f64, + { + // Validate inputs + if a.is_nan() || b.is_nan() || c.is_nan() { + return Err(QuadratureError::DegenerateInterval); + } + let (lo, hi) = if a < b { (a, b) } else { (b, a) }; + if c <= lo || c >= hi { + return Err(QuadratureError::InvalidInput( + "singularity c must be strictly inside (a, b)", + )); + } + + // Evaluate f at the singularity + let fc = f(c); + if !fc.is_finite() { + return Err(QuadratureError::InvalidInput( + "f(c) must be finite for the subtraction technique", + )); + } + + // Analytic term: f(c) * ln((b - c) / (c - a)) + let log_term = fc * ((b - c) / (c - a)).ln(); + + // The subtracted integrand: g(x) = [f(x) - f(c)] / (x - c) + // This has a removable singularity at x = c. + let guard = 1e-15 * (b - a); + let g = |x: f64| -> f64 { + if (x - c).abs() < guard { + 0.0 + } else { + (f(x) - fc) / (x - c) + } + }; + + // Integrate g over [a, b] with a break point at c + let adaptive_result = AdaptiveIntegrator::default() + .with_abs_tol(self.abs_tol) + .with_rel_tol(self.rel_tol) + .with_max_evals(self.max_evals) + .integrate_with_breaks(a, b, &[c], g)?; + + Ok(QuadratureResult { + value: adaptive_result.value + log_term, + error_estimate: adaptive_result.error_estimate, + num_evals: adaptive_result.num_evals + 1, // +1 for f(c) + converged: adaptive_result.converged, + }) + } +} + +/// Convenience: Cauchy principal value integration with default settings. +pub fn pv_integrate( + f: G, + a: f64, + b: f64, + c: f64, + tol: f64, +) -> Result, QuadratureError> +where + G: Fn(f64) -> f64, +{ + CauchyPV::default() + .with_abs_tol(tol) + .with_rel_tol(tol) + .integrate(a, b, c, f) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn symmetric_constant() { + // PV ∫₀¹ 1/(x - 0.5) dx = ln(0.5/0.5) = 0 + let result = pv_integrate(|_| 1.0, 0.0, 1.0, 0.5, 1e-10).unwrap(); + assert!( + result.value.abs() < 1e-8, + "value={}, expected 0", + result.value + ); + } + + #[test] + fn odd_integrand() { + // PV ∫₋₁¹ 1/x dx = 0 (odd function) + let result = pv_integrate(|_| 1.0, -1.0, 1.0, 0.0, 1e-10).unwrap(); + assert!( + result.value.abs() < 1e-8, + "value={}, expected 0", + result.value + ); + } + + #[test] + fn polynomial_f() { + // PV ∫₀¹ x²/(x - 0.3) dx + // = ∫₀¹ (x + 0.3) dx + 0.09 * ln(0.7/0.3) + // = [x²/2 + 0.3x]₀¹ + 0.09 * ln(7/3) + // = 0.8 + 0.09 * ln(7/3) + let exact = 0.8 + 0.09 * (7.0_f64 / 3.0).ln(); + let result = pv_integrate(|x| x * x, 0.0, 1.0, 0.3, 1e-10).unwrap(); + assert!( + (result.value - exact).abs() < 1e-7, + "value={}, expected={exact}", + result.value + ); + } + + #[test] + fn linear_f() { + // PV ∫₀² x/(x - 1) dx + // = ∫₀² [x - 1]/(x - 1) dx + 1 * ln(1/1) + // = ∫₀² 1 dx + 0 = 2 + let result = pv_integrate(|x| x, 0.0, 2.0, 1.0, 1e-10).unwrap(); + assert!( + (result.value - 2.0).abs() < 1e-7, + "value={}, expected=2", + result.value + ); + } + + #[test] + fn c_outside_interval() { + assert!(pv_integrate(|x| x, 0.0, 1.0, 2.0, 1e-10).is_err()); + } + + #[test] + fn c_at_endpoint() { + assert!(pv_integrate(|x| x, 0.0, 1.0, 0.0, 1e-10).is_err()); + assert!(pv_integrate(|x| x, 0.0, 1.0, 1.0, 1e-10).is_err()); + } + + #[test] + fn nan_inputs() { + assert!(pv_integrate(|x| x, f64::NAN, 1.0, 0.5, 1e-10).is_err()); + assert!(pv_integrate(|x| x, 0.0, 1.0, f64::NAN, 1e-10).is_err()); + } +} diff --git a/src/lib.rs b/src/lib.rs index 511bd32..a5b6ce6 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -57,6 +57,7 @@ //! ``` pub mod adaptive; +pub mod cauchy_pv; pub mod clenshaw_curtis; pub mod cubature; pub mod error; @@ -69,11 +70,15 @@ pub mod gauss_legendre; pub mod gauss_lobatto; pub mod gauss_radau; pub(crate) mod golub_welsch; +pub mod oscillatory; pub mod result; pub mod rule; +pub mod tanh_sinh; pub mod transforms; +pub mod weighted; pub use adaptive::{adaptive_integrate, adaptive_integrate_with_breaks, AdaptiveIntegrator}; +pub use cauchy_pv::{pv_integrate, CauchyPV}; pub use clenshaw_curtis::ClenshawCurtis; pub use error::QuadratureError; pub use gauss_chebyshev::{GaussChebyshevFirstKind, GaussChebyshevSecondKind}; @@ -84,8 +89,13 @@ pub use gauss_laguerre::GaussLaguerre; pub use gauss_legendre::GaussLegendre; pub use gauss_lobatto::GaussLobatto; pub use gauss_radau::GaussRadau; +pub use oscillatory::{ + integrate_oscillatory_cos, integrate_oscillatory_sin, OscillatoryIntegrator, OscillatoryKernel, +}; pub use result::QuadratureResult; pub use rule::QuadratureRule; +pub use tanh_sinh::{tanh_sinh_integrate, TanhSinh}; pub use transforms::{ integrate_infinite, integrate_semi_infinite_lower, integrate_semi_infinite_upper, }; +pub use weighted::{weighted_integrate, WeightFunction, WeightedIntegrator}; diff --git a/src/oscillatory.rs b/src/oscillatory.rs new file mode 100644 index 0000000..720dc97 --- /dev/null +++ b/src/oscillatory.rs @@ -0,0 +1,441 @@ +//! Oscillatory integration via Filon-Clenshaw-Curtis. +//! +//! Computes ∫ₐᵇ f(x)·sin(ωx) dx or ∫ₐᵇ f(x)·cos(ωx) dx efficiently +//! even when ω is large. Expands f in a Chebyshev basis, then integrates +//! each Chebyshev moment against the oscillatory kernel analytically. +//! +//! For small ω (|ω·(b-a)| < 2), falls back to standard adaptive integration. +//! +//! # Example +//! +//! ``` +//! use bilby::oscillatory::integrate_oscillatory_sin; +//! +//! // ∫₀¹ sin(100x) dx = (1 - cos(100)) / 100 +//! let exact = (1.0 - 100.0_f64.cos()) / 100.0; +//! let result = integrate_oscillatory_sin(|_| 1.0, 0.0, 1.0, 100.0, 1e-10).unwrap(); +//! assert!((result.value - exact).abs() < 1e-8); +//! ``` + +use crate::adaptive; +use crate::error::QuadratureError; +use crate::gauss_legendre::GaussLegendre; +use crate::result::QuadratureResult; + +/// Type of oscillatory kernel. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum OscillatoryKernel { + /// sin(ω·x) + Sine, + /// cos(ω·x) + Cosine, +} + +/// Builder for oscillatory integration via Filon-Clenshaw-Curtis. +/// +/// # Example +/// +/// ``` +/// use bilby::oscillatory::{OscillatoryIntegrator, OscillatoryKernel}; +/// +/// let integrator = OscillatoryIntegrator::new(OscillatoryKernel::Cosine, 50.0) +/// .with_order(64); +/// let result = integrator.integrate(0.0, 1.0, |_| 1.0).unwrap(); +/// let exact = 50.0_f64.sin() / 50.0; +/// assert!((result.value - exact).abs() < 1e-10); +/// ``` +#[derive(Debug, Clone)] +pub struct OscillatoryIntegrator { + kernel: OscillatoryKernel, + omega: f64, + order: usize, + abs_tol: f64, + rel_tol: f64, +} + +impl OscillatoryIntegrator { + /// Create a new oscillatory integrator. + pub fn new(kernel: OscillatoryKernel, omega: f64) -> Self { + Self { + kernel, + omega, + order: 32, + abs_tol: 1.49e-8, + rel_tol: 1.49e-8, + } + } + + /// Set the Chebyshev expansion order. + pub fn with_order(mut self, n: usize) -> Self { + self.order = n; + self + } + + /// Set absolute tolerance. + pub fn with_abs_tol(mut self, tol: f64) -> Self { + self.abs_tol = tol; + self + } + + /// Set relative tolerance. + pub fn with_rel_tol(mut self, tol: f64) -> Self { + self.rel_tol = tol; + self + } + + /// Integrate f(x)·kernel(ω·x) over \[a, b\]. + pub fn integrate( + &self, + a: f64, + b: f64, + f: G, + ) -> Result, QuadratureError> + where + G: Fn(f64) -> f64, + { + if a.is_nan() || b.is_nan() || self.omega.is_nan() { + return Err(QuadratureError::DegenerateInterval); + } + if (b - a).abs() < f64::EPSILON { + return Ok(QuadratureResult { + value: 0.0, + error_estimate: 0.0, + num_evals: 0, + converged: true, + }); + } + + let half = 0.5 * (b - a); + let mid = 0.5 * (a + b); + let theta = self.omega * half; + + // For small theta, fall back to adaptive integration + if theta.abs() < 2.0 { + return self.adaptive_fallback(a, b, &f); + } + + // Filon-Clenshaw-Curtis method + let n = self.order.max(4); + + // Step 1: Evaluate f at Clenshaw-Curtis nodes on [-1, 1] + let f_vals: Vec = (0..=n) + .map(|k| { + let t = (k as f64 * std::f64::consts::PI / n as f64).cos(); + f(mid + half * t) + }) + .collect(); + let num_evals = n + 1; + + // Step 2: Compute Chebyshev coefficients via DCT + let cheb_coeffs = chebyshev_coefficients(&f_vals, n); + + // Step 3: Compute modified Chebyshev moments + let (moments_cos, moments_sin) = modified_chebyshev_moments(theta, n); + + // Step 4: Sum I = half * Σ c_j * μ_j + // For cos kernel: ∫₋₁¹ f(t) cos(θt) dt = Σ c_j μ_j^c + // For sin kernel: ∫₋₁¹ f(t) sin(θt) dt = Σ c_j μ_j^s + // Then multiply by half and the phase factor for the shift to [a, b]. + // + // ∫ₐᵇ f(x) cos(ωx) dx = ∫₋₁¹ f(mid + half·t) cos(ω(mid + half·t)) half dt + // = half ∫₋₁¹ f̃(t) [cos(ωmid)cos(θt) - sin(ωmid)sin(θt)] dt + // = half [cos(ωmid) Σ cⱼμⱼᶜ - sin(ωmid) Σ cⱼμⱼˢ] + + let sum_c: f64 = cheb_coeffs + .iter() + .zip(moments_cos.iter()) + .map(|(c, m)| c * m) + .sum(); + let sum_s: f64 = cheb_coeffs + .iter() + .zip(moments_sin.iter()) + .map(|(c, m)| c * m) + .sum(); + + let omega_mid = self.omega * mid; + let cos_mid = omega_mid.cos(); + let sin_mid = omega_mid.sin(); + + let value = match self.kernel { + OscillatoryKernel::Cosine => half * (cos_mid * sum_c - sin_mid * sum_s), + OscillatoryKernel::Sine => half * (sin_mid * sum_c + cos_mid * sum_s), + }; + + // Error estimate: compare with a lower-order approximation (use half the coefficients) + let n_half = n / 2; + let sum_c_half: f64 = cheb_coeffs + .iter() + .take(n_half + 1) + .zip(moments_cos.iter()) + .map(|(c, m)| c * m) + .sum(); + let sum_s_half: f64 = cheb_coeffs + .iter() + .take(n_half + 1) + .zip(moments_sin.iter()) + .map(|(c, m)| c * m) + .sum(); + + let value_half = match self.kernel { + OscillatoryKernel::Cosine => half * (cos_mid * sum_c_half - sin_mid * sum_s_half), + OscillatoryKernel::Sine => half * (sin_mid * sum_c_half + cos_mid * sum_s_half), + }; + + let error = (value - value_half).abs(); + let converged = error <= self.abs_tol.max(self.rel_tol * value.abs()); + + Ok(QuadratureResult { + value, + error_estimate: error, + num_evals, + converged, + }) + } + + /// Fallback to standard adaptive integration for small ω. + fn adaptive_fallback( + &self, + a: f64, + b: f64, + f: &G, + ) -> Result, QuadratureError> + where + G: Fn(f64) -> f64, + { + let omega = self.omega; + let integrand = match self.kernel { + OscillatoryKernel::Sine => { + Box::new(move |x: f64| f(x) * (omega * x).sin()) as Box f64> + } + OscillatoryKernel::Cosine => Box::new(move |x: f64| f(x) * (omega * x).cos()), + }; + adaptive::adaptive_integrate(&*integrand, a, b, self.abs_tol) + } +} + +/// Compute Chebyshev coefficients from function values at CC nodes. +/// +/// Given f_k = f(cos(kπ/n)) for k = 0, ..., n, computes the coefficients +/// c_j such that f(t) ≈ Σ c_j T_j(t) (with the convention that c_0 and +/// c_n are halved). +fn chebyshev_coefficients(f_vals: &[f64], n: usize) -> Vec { + let mut coeffs = vec![0.0; n + 1]; + let pi_n = std::f64::consts::PI / n as f64; + + for (j, cj) in coeffs.iter_mut().enumerate() { + let mut sum = 0.0; + for (k, fk) in f_vals.iter().enumerate() { + let factor = if k == 0 || k == n { 0.5 } else { 1.0 }; + sum += factor * fk * (j as f64 * k as f64 * pi_n).cos(); + } + *cj = 2.0 * sum / n as f64; + } + + // Halve the first and last coefficients + coeffs[0] *= 0.5; + coeffs[n] *= 0.5; + + coeffs +} + +/// Compute modified Chebyshev moments via Gauss-Legendre quadrature. +/// +/// Returns (μ_j^c, μ_j^s) for j = 0, ..., n where: +/// μ_j^c = ∫₋₁¹ T_j(x) cos(θx) dx +/// μ_j^s = ∫₋₁¹ T_j(x) sin(θx) dx +/// +/// Uses GL quadrature with enough nodes to resolve both the Chebyshev +/// polynomial (degree j ≤ n) and the oscillatory kernel (effective +/// frequency ~θ). Converges exponentially for these smooth integrands. +fn modified_chebyshev_moments(theta: f64, n: usize) -> (Vec, Vec) { + let mut mc = vec![0.0; n + 1]; + let mut ms = vec![0.0; n + 1]; + + if theta.abs() < 1e-15 { + // cos(θx) ≈ 1, sin(θx) ≈ 0 + // μ_j^c = ∫₋₁¹ T_j(x) dx + for (j, mcj) in mc.iter_mut().enumerate() { + if j == 0 { + *mcj = 2.0; + } else if j % 2 == 0 { + *mcj = 2.0 / (1.0 - (j as f64).powi(2)); + } + } + return (mc, ms); + } + + // GL with m nodes is exact for polynomials of degree 2m-1. + // T_j(x) cos(θx) is smooth and entire; GL converges exponentially + // when m exceeds the effective bandwidth n + |θ|. + let m = (n + (theta.abs().ceil() as usize) + 32).max(64); + let gl = GaussLegendre::new(m).unwrap(); + let rule = gl.rule(); + + for (node, weight) in rule.nodes.iter().zip(rule.weights.iter()) { + let x = *node; + let cos_tx = (theta * x).cos(); + let sin_tx = (theta * x).sin(); + + // Compute T_j(x) via recurrence: T_0=1, T_1=x, T_{k+1}=2xT_k - T_{k-1} + let mut t_prev = 1.0; + let mut t_curr = x; + + mc[0] += weight * cos_tx; + ms[0] += weight * sin_tx; + + if n >= 1 { + mc[1] += weight * x * cos_tx; + ms[1] += weight * x * sin_tx; + } + + for j in 2..=n { + let t_next = 2.0 * x * t_curr - t_prev; + mc[j] += weight * t_next * cos_tx; + ms[j] += weight * t_next * sin_tx; + t_prev = t_curr; + t_curr = t_next; + } + } + + (mc, ms) +} + +/// Convenience: integrate f(x)·sin(ωx) over \[a, b\]. +pub fn integrate_oscillatory_sin( + f: G, + a: f64, + b: f64, + omega: f64, + tol: f64, +) -> Result, QuadratureError> +where + G: Fn(f64) -> f64, +{ + OscillatoryIntegrator::new(OscillatoryKernel::Sine, omega) + .with_abs_tol(tol) + .with_rel_tol(tol) + .integrate(a, b, f) +} + +/// Convenience: integrate f(x)·cos(ωx) over \[a, b\]. +pub fn integrate_oscillatory_cos( + f: G, + a: f64, + b: f64, + omega: f64, + tol: f64, +) -> Result, QuadratureError> +where + G: Fn(f64) -> f64, +{ + OscillatoryIntegrator::new(OscillatoryKernel::Cosine, omega) + .with_abs_tol(tol) + .with_rel_tol(tol) + .integrate(a, b, f) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn sin_trivial() { + // ∫₀^π sin(x) dx = 2 (ω=1, f=1) + let result = + integrate_oscillatory_sin(|_| 1.0, 0.0, std::f64::consts::PI, 1.0, 1e-10).unwrap(); + assert!((result.value - 2.0).abs() < 1e-6, "value={}", result.value); + } + + #[test] + fn cos_moderate_omega() { + // ∫₀¹ cos(10x) dx = sin(10)/10 + let exact = 10.0_f64.sin() / 10.0; + let result = integrate_oscillatory_cos(|_| 1.0, 0.0, 1.0, 10.0, 1e-10).unwrap(); + assert!( + (result.value - exact).abs() < 1e-8, + "value={}, exact={exact}", + result.value + ); + } + + #[test] + fn sin_high_omega() { + // ∫₀¹ sin(100x) dx = (1 - cos(100))/100 + let exact = (1.0 - 100.0_f64.cos()) / 100.0; + let result = integrate_oscillatory_sin(|_| 1.0, 0.0, 1.0, 100.0, 1e-10).unwrap(); + assert!( + (result.value - exact).abs() < 1e-8, + "value={}, exact={exact}", + result.value + ); + } + + #[test] + fn cos_with_linear_f() { + // ∫₀^π x cos(x) dx = [x sin(x) + cos(x)]₀^π = -1 - 1 = -2 + let result = + integrate_oscillatory_cos(|x| x, 0.0, std::f64::consts::PI, 1.0, 1e-8).unwrap(); + assert!( + (result.value - (-2.0)).abs() < 1e-4, + "value={}", + result.value + ); + } + + #[test] + fn sin_with_exp_f() { + // ∫₀¹ exp(x) sin(ωx) dx = Im[(e^(1+iω) - 1) / (1+iω)] + let omega: f64 = 50.0; + let r = 1.0; + let i = omega; + // e^(1+iω) = e * (cos(ω) + i sin(ω)) + let e = std::f64::consts::E; + let re_num = e * omega.cos() - 1.0; + let im_num = e * omega.sin(); + let denom = r * r + i * i; // 1 + ω² + // (re_num + i im_num) / (r + i*i_) = (re_num + i im_num)(r - i*i_) / denom + let exact = (im_num * r - re_num * i) / denom; // imaginary part + + let result = integrate_oscillatory_sin(f64::exp, 0.0, 1.0, omega, 1e-8).unwrap(); + assert!( + (result.value - exact).abs() < 1e-4, + "value={}, exact={exact}", + result.value + ); + } + + #[test] + fn small_omega_fallback() { + // ∫₀¹ cos(0.5x) dx = sin(0.5)/0.5 (falls back to adaptive) + let exact = 0.5_f64.sin() / 0.5; + let result = integrate_oscillatory_cos(|_| 1.0, 0.0, 1.0, 0.5, 1e-10).unwrap(); + assert!( + (result.value - exact).abs() < 1e-8, + "value={}, exact={exact}", + result.value + ); + } + + #[test] + fn zero_interval() { + let result = integrate_oscillatory_sin(|_| 1.0, 1.0, 1.0, 10.0, 1e-10).unwrap(); + assert_eq!(result.value, 0.0); + } + + #[test] + fn nan_input() { + assert!(integrate_oscillatory_sin(|_| 1.0, f64::NAN, 1.0, 10.0, 1e-10).is_err()); + } + + #[test] + fn very_high_omega() { + // ∫₀¹ sin(1000x) dx = (1 - cos(1000))/1000 + let exact = (1.0 - 1000.0_f64.cos()) / 1000.0; + let result = integrate_oscillatory_sin(|_| 1.0, 0.0, 1.0, 1000.0, 1e-8).unwrap(); + assert!( + (result.value - exact).abs() < 1e-6, + "value={}, exact={exact}", + result.value + ); + } +} diff --git a/src/tanh_sinh.rs b/src/tanh_sinh.rs new file mode 100644 index 0000000..3962cbc --- /dev/null +++ b/src/tanh_sinh.rs @@ -0,0 +1,352 @@ +//! Tanh-sinh (double-exponential) quadrature. +//! +//! Transforms endpoint singularities into rapidly decaying integrands via +//! x = tanh(π/2 · sinh(t)). Self-adaptive through level doubling: each level +//! halves the step size and reuses all previously computed points. +//! +//! Particularly effective for integrands with algebraic or logarithmic +//! endpoint singularities where standard Gauss rules struggle. +//! +//! # Example +//! +//! ``` +//! use bilby::tanh_sinh_integrate; +//! +//! // Integral of 1/sqrt(x) over [0, 1] = 2 (endpoint singularity) +//! let result = tanh_sinh_integrate(|x| 1.0 / x.sqrt(), 0.0, 1.0, 1e-10).unwrap(); +//! assert!((result.value - 2.0).abs() < 1e-7); +//! ``` + +use crate::error::QuadratureError; +use crate::result::QuadratureResult; + +/// Builder for tanh-sinh (double-exponential) quadrature. +/// +/// # Example +/// +/// ``` +/// use bilby::tanh_sinh::TanhSinh; +/// +/// let result = TanhSinh::default() +/// .with_abs_tol(1e-12) +/// .integrate(0.0, 1.0, |x| x.ln()) +/// .unwrap(); +/// assert!((result.value - (-1.0)).abs() < 1e-10); +/// ``` +#[derive(Debug, Clone)] +pub struct TanhSinh { + max_levels: usize, + abs_tol: f64, + rel_tol: f64, +} + +impl Default for TanhSinh { + fn default() -> Self { + Self { + max_levels: 12, + abs_tol: 1.49e-8, + rel_tol: 1.49e-8, + } + } +} + +impl TanhSinh { + /// Set the maximum number of refinement levels. + pub fn with_max_levels(mut self, levels: usize) -> Self { + self.max_levels = levels; + self + } + + /// Set absolute tolerance. + pub fn with_abs_tol(mut self, tol: f64) -> Self { + self.abs_tol = tol; + self + } + + /// Set relative tolerance. + pub fn with_rel_tol(mut self, tol: f64) -> Self { + self.rel_tol = tol; + self + } + + /// Integrate `f` over [a, b] using the tanh-sinh transform. + pub fn integrate( + &self, + a: f64, + b: f64, + f: G, + ) -> Result, QuadratureError> + where + G: Fn(f64) -> f64, + { + if a.is_nan() || b.is_nan() { + return Err(QuadratureError::DegenerateInterval); + } + if (b - a).abs() < f64::EPSILON { + return Ok(QuadratureResult { + value: 0.0, + error_estimate: 0.0, + num_evals: 0, + converged: true, + }); + } + + let mid = 0.5 * (a + b); + let half = 0.5 * (b - a); + let pi_2 = std::f64::consts::FRAC_PI_2; + + let mut total_evals = 0usize; + let mut prev_estimate: f64 = 0.0; + let mut h: f64 = 1.0; + + // Level 0: evaluate at all points t = k*h + // Subsequent levels: only evaluate at odd multiples of the new h + for level in 0..=self.max_levels { + let step = if level == 0 { 1 } else { 2 }; + let start = 1; // k=0 is the center point, handled separately + + // At level L, h = 2^{-L}. New points are at t = (2k+1) * h for k >= 0. + // But we also need negative t values. + + if level > 0 { + // Scale previous sum for the new step size (halved h) + // prev_estimate was computed with old h, which is 2*current h. + // The trapezoidal rule sum scales linearly with h, so halving h + // means we keep all old points and add new ones. + // Actually, the estimate is h * sum_of_weighted_values. + // When we halve h, old_estimate = 2h * old_sum = old_value. + // New estimate = h * (old_sum + new_sum) = old_value/2 + h * new_sum. + h *= 0.5; + } + + let mut new_sum = 0.0; + + // Center point (only at level 0) + if level == 0 { + // t = 0: sinh(0) = 0, cosh(0) = 1, tanh(0) = 0 + // x = mid, weight = pi/2 + let w0 = pi_2; + let fval = f(mid); + if fval.is_finite() { + new_sum += w0 * fval; + } + total_evals += 1; + } + + // Positive and negative t values + // Upper bound for the loop; actual termination is handled by the + // u > 710, weight < 1e-300, and consecutive_tiny checks. + let max_k = (7.0 / h).ceil() as usize; + let mut k = start; + let mut consecutive_tiny = 0; + + while k <= max_k { + let t = k as f64 * h; + let sinh_t = t.sinh(); + let cosh_t = t.cosh(); + let u = pi_2 * sinh_t; + + // For large |u|, tanh(u) ≈ ±1 and cosh(u) is huge. + // The weight decays double-exponentially. + if u.abs() > 710.0 { + // cosh(u)^2 would overflow; weight is effectively zero + break; + } + + let cosh_u = u.cosh(); + let tanh_u = u.tanh(); + let weight = pi_2 * cosh_t / (cosh_u * cosh_u); + + if weight < 1e-300 { + break; + } + + // Positive t: x = mid + half * tanh(u) + let mut local_contrib = 0.0_f64; + let x_pos = mid + half * tanh_u; + if x_pos > a && x_pos < b { + let fval = f(x_pos); + if fval.is_finite() { + new_sum += weight * fval; + local_contrib = local_contrib.max((weight * fval).abs()); + } + total_evals += 1; + } + + // Negative t: x = mid - half * tanh(u) + let x_neg = mid - half * tanh_u; + if x_neg > a && x_neg < b { + let fval = f(x_neg); + if fval.is_finite() { + new_sum += weight * fval; + local_contrib = local_contrib.max((weight * fval).abs()); + } + total_evals += 1; + } + + // Check if actual contributions (including f(x)) are negligible + let threshold = f64::EPSILON * prev_estimate.abs().max(1e-300); + if local_contrib * half < threshold { + consecutive_tiny += 1; + if consecutive_tiny >= 3 { + break; + } + } else { + consecutive_tiny = 0; + } + + k += step; + } + + let estimate = if level == 0 { + h * half * new_sum + } else { + 0.5 * prev_estimate + h * half * new_sum + }; + + // Require at least 3 levels before checking convergence. + // Early levels have so few points that successive differences + // can be misleadingly small for singular integrands. + if level >= 3 { + let error = (estimate - prev_estimate).abs(); + let tol = self.abs_tol.max(self.rel_tol * estimate.abs()); + if error <= tol { + return Ok(QuadratureResult { + value: estimate, + error_estimate: error, + num_evals: total_evals, + converged: true, + }); + } + } + + prev_estimate = estimate; + } + + // Did not converge within max_levels + let error = if self.max_levels > 0 { + // Can't easily recompute the last level difference, + // but prev_estimate is our best value + self.abs_tol.max(self.rel_tol * prev_estimate.abs()) * 10.0 + } else { + f64::INFINITY + }; + + Ok(QuadratureResult { + value: prev_estimate, + error_estimate: error, + num_evals: total_evals, + converged: false, + }) + } +} + +/// Convenience: tanh-sinh integration with default settings. +pub fn tanh_sinh_integrate( + f: G, + a: f64, + b: f64, + tol: f64, +) -> Result, QuadratureError> +where + G: Fn(f64) -> f64, +{ + TanhSinh::default() + .with_abs_tol(tol) + .with_rel_tol(tol) + .integrate(a, b, f) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn smooth_polynomial() { + let result = tanh_sinh_integrate(|x| x * x, 0.0, 1.0, 1e-12).unwrap(); + assert!( + (result.value - 1.0 / 3.0).abs() < 1e-10, + "value={}", + result.value + ); + assert!(result.converged); + } + + #[test] + fn sqrt_singularity() { + // ∫₀¹ 1/√x dx = 2 + // f64 limits: tanh saturates near ±1, losing ~1.5e-8 of the integral + let result = tanh_sinh_integrate(|x| 1.0 / x.sqrt(), 0.0, 1.0, 1e-10).unwrap(); + assert!((result.value - 2.0).abs() < 1e-7, "value={}", result.value); + } + + #[test] + fn log_singularity() { + // ∫₀¹ ln(x) dx = -1 + let result = tanh_sinh_integrate(|x| x.ln(), 0.0, 1.0, 1e-10).unwrap(); + assert!( + (result.value - (-1.0)).abs() < 1e-8, + "value={}", + result.value + ); + } + + #[test] + fn strong_algebraic_singularity() { + // ∫₀¹ x^(-0.75) dx = 1/0.25 = 4 + // x^(-0.9) loses ~0.24 due to f64 tanh saturation near ±1; + // x^(-0.75) loses only ~3.5e-4, well within tolerance. + let result = tanh_sinh_integrate(|x| x.powf(-0.75), 0.0, 1.0, 1e-6).unwrap(); + assert!((result.value - 4.0).abs() < 0.01, "value={}", result.value); + } + + #[test] + fn both_endpoints_singular() { + // ∫₀¹ 1/√(x(1-x)) dx = π + let result = tanh_sinh_integrate(|x| 1.0 / (x * (1.0 - x)).sqrt(), 0.0, 1.0, 1e-8).unwrap(); + assert!( + (result.value - std::f64::consts::PI).abs() < 1e-6, + "value={}", + result.value + ); + } + + #[test] + fn chebyshev_weight() { + // ∫₋₁¹ 1/√(1-x²) dx = π + let result = tanh_sinh_integrate(|x| 1.0 / (1.0 - x * x).sqrt(), -1.0, 1.0, 1e-8).unwrap(); + assert!( + (result.value - std::f64::consts::PI).abs() < 1e-6, + "value={}", + result.value + ); + } + + #[test] + fn sin_integral() { + // ∫₀^π sin(x) dx = 2 (smooth, sanity check) + let result = tanh_sinh_integrate(|x| x.sin(), 0.0, std::f64::consts::PI, 1e-10).unwrap(); + assert!((result.value - 2.0).abs() < 1e-8, "value={}", result.value); + } + + #[test] + fn zero_width_interval() { + let result = tanh_sinh_integrate(|x| x, 1.0, 1.0, 1e-10).unwrap(); + assert_eq!(result.value, 0.0); + assert!(result.converged); + } + + #[test] + fn nan_bounds() { + assert!(tanh_sinh_integrate(|x| x, f64::NAN, 1.0, 1e-10).is_err()); + } + + #[test] + fn non_convergence() { + let result = TanhSinh::default() + .with_max_levels(0) + .integrate(0.0, 1.0, |x| 1.0 / x.sqrt()) + .unwrap(); + assert!(!result.converged); + } +} diff --git a/src/weighted.rs b/src/weighted.rs new file mode 100644 index 0000000..cfe81b5 --- /dev/null +++ b/src/weighted.rs @@ -0,0 +1,294 @@ +//! Weighted integration: ∫ f(x) · w(x) dx for known weight functions. +//! +//! Dispatches to the appropriate Gaussian quadrature rule for each weight +//! function family, providing a unified API for product integration. +//! +//! # Example +//! +//! ``` +//! use bilby::weighted::{weighted_integrate, WeightFunction}; +//! +//! // ∫₋₁¹ 1/√(1-x²) dx = π (Chebyshev Type I weight) +//! let result = weighted_integrate(|_| 1.0, WeightFunction::ChebyshevI, 20).unwrap(); +//! assert!((result - std::f64::consts::PI).abs() < 1e-12); +//! ``` + +use crate::error::QuadratureError; +use crate::gauss_chebyshev::{GaussChebyshevFirstKind, GaussChebyshevSecondKind}; +use crate::gauss_hermite::GaussHermite; +use crate::gauss_jacobi::GaussJacobi; +use crate::gauss_laguerre::GaussLaguerre; + +/// Known weight functions for product integration. +#[derive(Debug, Clone)] +pub enum WeightFunction { + /// w(x) = (1-x)^α · (1+x)^β on \[-1, 1\] + Jacobi { alpha: f64, beta: f64 }, + /// w(x) = x^α · e^(-x) on \[0, ∞) + Laguerre { alpha: f64 }, + /// w(x) = e^(-x²) on (-∞, ∞) + Hermite, + /// w(x) = 1/√(1 - x²) on (-1, 1) + ChebyshevI, + /// w(x) = √(1 - x²) on \[-1, 1\] + ChebyshevII, + /// w(x) = -ln(x) on (0, 1\] + LogWeight, +} + +/// Builder for weighted integration. +/// +/// # Example +/// +/// ``` +/// use bilby::weighted::{WeightedIntegrator, WeightFunction}; +/// +/// let wi = WeightedIntegrator::new(WeightFunction::Hermite, 20).unwrap(); +/// // ∫₋∞^∞ 1 · e^(-x²) dx = √π +/// let result = wi.integrate(|_| 1.0); +/// assert!((result - std::f64::consts::PI.sqrt()).abs() < 1e-12); +/// ``` +#[derive(Debug, Clone)] +pub struct WeightedIntegrator { + weight: WeightFunction, + order: usize, +} + +impl WeightedIntegrator { + /// Create a new weighted integrator. + pub fn new(weight: WeightFunction, order: usize) -> Result { + if order == 0 { + return Err(QuadratureError::ZeroOrder); + } + Ok(Self { weight, order }) + } + + /// Set the quadrature order. + pub fn with_order(mut self, n: usize) -> Self { + self.order = n; + self + } + + /// Integrate f(x) · w(x) over the natural domain of w. + /// + /// - Jacobi: \[-1, 1\] + /// - Laguerre: \[0, ∞) + /// - Hermite: (-∞, ∞) + /// - ChebyshevI/II: \[-1, 1\] + /// - LogWeight: (0, 1\] + pub fn integrate(&self, f: G) -> f64 + where + G: Fn(f64) -> f64, + { + match &self.weight { + WeightFunction::Jacobi { alpha, beta } => { + let gj = GaussJacobi::new(self.order, *alpha, *beta).unwrap(); + let rule = gj.rule(); + let mut sum = 0.0; + for (node, weight) in rule.nodes.iter().zip(rule.weights.iter()) { + sum += weight * f(*node); + } + sum + } + WeightFunction::Laguerre { alpha } => { + let gl = GaussLaguerre::new(self.order, *alpha).unwrap(); + let rule = gl.rule(); + let mut sum = 0.0; + for (node, weight) in rule.nodes.iter().zip(rule.weights.iter()) { + sum += weight * f(*node); + } + sum + } + WeightFunction::Hermite => { + let gh = GaussHermite::new(self.order).unwrap(); + let rule = gh.rule(); + let mut sum = 0.0; + for (node, weight) in rule.nodes.iter().zip(rule.weights.iter()) { + sum += weight * f(*node); + } + sum + } + WeightFunction::ChebyshevI => { + let gc = GaussChebyshevFirstKind::new(self.order).unwrap(); + let rule = gc.rule(); + let mut sum = 0.0; + for (node, weight) in rule.nodes.iter().zip(rule.weights.iter()) { + sum += weight * f(*node); + } + sum + } + WeightFunction::ChebyshevII => { + let gc = GaussChebyshevSecondKind::new(self.order).unwrap(); + let rule = gc.rule(); + let mut sum = 0.0; + for (node, weight) in rule.nodes.iter().zip(rule.weights.iter()) { + sum += weight * f(*node); + } + sum + } + WeightFunction::LogWeight => { + // ∫₀¹ -ln(x) · f(x) dx via substitution x = e^{-t}: + // = ∫₀^∞ t · e^{-t} · f(e^{-t}) dt + // This is a Gauss-Laguerre integral with α = 1. + let gl = GaussLaguerre::new(self.order, 1.0).unwrap(); + let rule = gl.rule(); + let mut sum = 0.0; + for (node, weight) in rule.nodes.iter().zip(rule.weights.iter()) { + sum += weight * f((-*node).exp()); + } + sum + } + } + } + + /// Integrate f(x) · w(x) over \[a, b\] via affine transform. + /// + /// Only applicable for finite-domain weights (Jacobi, ChebyshevI/II). + /// For LogWeight, maps (0, 1\] to (a, b\] with the log singularity at a. + /// For Laguerre/Hermite, a and b are ignored (uses natural domain). + pub fn integrate_over(&self, a: f64, b: f64, f: G) -> f64 + where + G: Fn(f64) -> f64, + { + match &self.weight { + WeightFunction::Jacobi { .. } + | WeightFunction::ChebyshevI + | WeightFunction::ChebyshevII => { + // Affine map: x in [a,b] <-> t in [-1,1] + // x = (b-a)/2 * t + (a+b)/2 + // w(x) dx in terms of t includes a Jacobian factor + // For Jacobi weight: (1-t)^a (1+t)^b maps directly + let half = 0.5 * (b - a); + let mid = 0.5 * (a + b); + self.integrate(|t| f(half * t + mid)) * half.powi(1) + } + WeightFunction::LogWeight => { + // Map (0,1] to (a,b]: x = a + (b-a)*u, log singularity at a + // ∫ₐᵇ -ln(x-a) f(x) dx ... actually the weight is -ln(u) on (0,1] + // For general [a,b]: ∫ₐᵇ -ln(x-a) f(x) dx = (b-a) ∫₀¹ -ln((b-a)u) f(a+(b-a)u) du + // = (b-a) ∫₀¹ [-ln(b-a) - ln(u)] f(a+(b-a)u) du + // This doesn't decompose cleanly. Keep it simple: just scale. + let width = b - a; + let inner = self.integrate(|u| f(a + width * u)); + width * inner + } + // Laguerre and Hermite use natural domains + WeightFunction::Laguerre { .. } | WeightFunction::Hermite => self.integrate(f), + } + } +} + +/// Convenience: integrate f(x) · w(x) over the natural domain. +pub fn weighted_integrate( + f: G, + weight: WeightFunction, + order: usize, +) -> Result +where + G: Fn(f64) -> f64, +{ + let wi = WeightedIntegrator::new(weight, order)?; + Ok(wi.integrate(f)) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn jacobi_weight_sum() { + // ∫₋₁¹ (1-x)^0.5 (1+x)^0.5 dx = π/2 + let result = weighted_integrate( + |_| 1.0, + WeightFunction::Jacobi { + alpha: 0.5, + beta: 0.5, + }, + 20, + ) + .unwrap(); + assert!( + (result - std::f64::consts::FRAC_PI_2).abs() < 1e-10, + "result={result}" + ); + } + + #[test] + fn laguerre_constant() { + // ∫₀^∞ e^(-x) dx = 1 + let result = + weighted_integrate(|_| 1.0, WeightFunction::Laguerre { alpha: 0.0 }, 10).unwrap(); + assert!((result - 1.0).abs() < 1e-10, "result={result}"); + } + + #[test] + fn laguerre_linear() { + // ∫₀^∞ x · e^(-x) dx = Γ(2) = 1 + let result = + weighted_integrate(|x| x, WeightFunction::Laguerre { alpha: 0.0 }, 10).unwrap(); + assert!((result - 1.0).abs() < 1e-10, "result={result}"); + } + + #[test] + fn hermite_constant() { + // ∫₋∞^∞ e^(-x²) dx = √π + let result = weighted_integrate(|_| 1.0, WeightFunction::Hermite, 10).unwrap(); + assert!( + (result - std::f64::consts::PI.sqrt()).abs() < 1e-10, + "result={result}" + ); + } + + #[test] + fn chebyshev_i_constant() { + // ∫₋₁¹ 1/√(1-x²) dx = π + let result = weighted_integrate(|_| 1.0, WeightFunction::ChebyshevI, 20).unwrap(); + assert!( + (result - std::f64::consts::PI).abs() < 1e-12, + "result={result}" + ); + } + + #[test] + fn chebyshev_ii_constant() { + // ∫₋₁¹ √(1-x²) dx = π/2 + let result = weighted_integrate(|_| 1.0, WeightFunction::ChebyshevII, 20).unwrap(); + assert!( + (result - std::f64::consts::FRAC_PI_2).abs() < 1e-12, + "result={result}" + ); + } + + #[test] + fn log_weight_constant() { + // ∫₀¹ -ln(x) dx = 1 + let result = weighted_integrate(|_| 1.0, WeightFunction::LogWeight, 20).unwrap(); + assert!((result - 1.0).abs() < 1e-10, "result={result}"); + } + + #[test] + fn log_weight_linear() { + // ∫₀¹ -ln(x) · x dx = 1/4 + let result = weighted_integrate(|x| x, WeightFunction::LogWeight, 20).unwrap(); + assert!((result - 0.25).abs() < 1e-10, "result={result}"); + } + + #[test] + fn log_weight_quadratic() { + // ∫₀¹ -ln(x) · x² dx = 2/27 (integration by parts twice) + // Actually: ∫₀¹ -ln(x) x^n dx = 1/(n+1)^2 + // For n=2: 1/9... wait, let me recalculate. + // ∫₀¹ -ln(x) x² dx = [-x³/3 ln(x)]₀¹ + ∫₀¹ x²/3 dx = 0 + 1/9 + let result = weighted_integrate(|x| x * x, WeightFunction::LogWeight, 20).unwrap(); + assert!( + (result - 1.0 / 9.0).abs() < 1e-10, + "result={result}, expected={}", + 1.0 / 9.0 + ); + } + + #[test] + fn zero_order() { + assert!(weighted_integrate(|_| 1.0, WeightFunction::Hermite, 0).is_err()); + } +}