bilby/
lib.rs

1#![cfg_attr(not(feature = "std"), no_std)]
2//! # bilby
3//!
4//! A high-performance numerical quadrature (integration) library for Rust.
5//!
6//! bilby provides Gaussian quadrature rules, adaptive integration, and
7//! multi-dimensional cubature methods with a consistent, generic API.
8//!
9//! ## Design
10//!
11//! - **Generic over `F: Float`** via [`num_traits::Float`] — works with `f32`, `f64`,
12//!   and compatible types (e.g., echidna AD types)
13//! - **Precomputed rules** — generate nodes and weights once, integrate many times
14//! - **No heap allocation on hot paths** — rule construction allocates, integration does not
15//! - **Separate 1D and N-D APIs** — `Fn(F) -> F` for 1D, `Fn(&[F]) -> F` for N-D
16//! - **`no_std` compatible** — works without the standard library (with `alloc`)
17//!
18//! ## Feature Flags
19//!
20//! | Feature | Default | Description |
21//! |---------|---------|-------------|
22//! | `std` | Yes | Enables `std::error::Error` impl and [`cache`] module |
23//! | `parallel` | No | Enables rayon-based `_par` methods (implies `std`) |
24//!
25//! ## Quick Start
26//!
27//! ```
28//! use bilby::GaussLegendre;
29//!
30//! // Create a 10-point Gauss-Legendre rule
31//! let gl = GaussLegendre::new(10).unwrap();
32//!
33//! // Integrate x^2 over [0, 1] (exact result = 1/3)
34//! let result = gl.rule().integrate(0.0, 1.0, |x: f64| x * x);
35//! assert!((result - 1.0 / 3.0).abs() < 1e-14);
36//! ```
37//!
38//! ## Gauss-Kronrod Error Estimation
39//!
40//! ```
41//! use bilby::{GaussKronrod, GKPair};
42//!
43//! let gk = GaussKronrod::new(GKPair::G7K15);
44//! let (estimate, error) = gk.integrate(0.0, core::f64::consts::PI, f64::sin);
45//! assert!((estimate - 2.0).abs() < 1e-14);
46//! ```
47//!
48//! ## Adaptive Integration
49//!
50//! ```
51//! use bilby::adaptive_integrate;
52//!
53//! let result = adaptive_integrate(|x: f64| x.sin(), 0.0, core::f64::consts::PI, 1e-12).unwrap();
54//! assert!((result.value - 2.0).abs() < 1e-12);
55//! assert!(result.is_converged());
56//! ```
57//!
58//! ## Infinite Domains
59//!
60//! ```
61//! use bilby::integrate_infinite;
62//!
63//! // Integral of exp(-x^2) over (-inf, inf) = sqrt(pi)
64//! let result = integrate_infinite(|x: f64| (-x * x).exp(), 1e-10).unwrap();
65//! assert!((result.value - core::f64::consts::PI.sqrt()).abs() < 1e-8);
66//! ```
67//!
68//! ## Multi-Dimensional Integration
69//!
70//! ```
71//! use bilby::cubature::{TensorProductRule, adaptive_cubature};
72//! use bilby::GaussLegendre;
73//!
74//! // Tensor product: 10-point GL in each of 2 dimensions
75//! let gl = GaussLegendre::new(10).unwrap();
76//! let tp = TensorProductRule::isotropic(gl.rule(), 2).unwrap();
77//! let result = tp.rule().integrate_box(
78//!     &[0.0, 0.0], &[1.0, 1.0],
79//!     |x| x[0] * x[1],
80//! );
81//! assert!((result - 0.25).abs() < 1e-14);
82//! ```
83//!
84//! ## Precomputed Rule Cache
85//!
86//! Available when the `std` feature is enabled (default):
87//!
88//! ```
89//! # #[cfg(feature = "std")] {
90//! use bilby::cache::GL10;
91//!
92//! let result = GL10.rule().integrate(0.0, 1.0, |x: f64| x * x);
93//! assert!((result - 1.0 / 3.0).abs() < 1e-14);
94//! # }
95//! ```
96//!
97//! ## Tanh-Sinh (Double Exponential)
98//!
99//! ```
100//! use bilby::tanh_sinh_integrate;
101//!
102//! // Endpoint singularity: integral of 1/sqrt(x) over [0, 1] = 2
103//! let result = tanh_sinh_integrate(|x| 1.0 / x.sqrt(), 0.0, 1.0, 1e-10).unwrap();
104//! assert!((result.value - 2.0).abs() < 1e-7);
105//! ```
106//!
107//! ## Cauchy Principal Value
108//!
109//! ```
110//! use bilby::pv_integrate;
111//!
112//! // PV integral of x^2/(x - 0.3) over [0, 1]
113//! let exact = 0.8 + 0.09 * (7.0_f64 / 3.0).ln();
114//! let result = pv_integrate(|x| x * x, 0.0, 1.0, 0.3, 1e-10).unwrap();
115//! assert!((result.value - exact).abs() < 1e-7);
116//! ```
117//!
118//! ## Oscillatory Integration
119//!
120//! ```
121//! use bilby::integrate_oscillatory_sin;
122//!
123//! // Integral of sin(100x) over [0, 1]
124//! let exact = (1.0 - 100.0_f64.cos()) / 100.0;
125//! let result = integrate_oscillatory_sin(|_| 1.0, 0.0, 1.0, 100.0, 1e-10).unwrap();
126//! assert!((result.value - exact).abs() < 1e-8);
127//! ```
128//!
129//! ## Weighted Integration
130//!
131//! ```
132//! use bilby::weighted::{weighted_integrate, WeightFunction};
133//!
134//! // Gauss-Hermite: integral of e^(-x^2) over (-inf, inf) = sqrt(pi)
135//! let result = weighted_integrate(|_| 1.0, WeightFunction::Hermite, 20).unwrap();
136//! assert!((result - core::f64::consts::PI.sqrt()).abs() < 1e-12);
137//! ```
138//!
139//! ## Sparse Grids
140//!
141//! ```
142//! use bilby::cubature::SparseGrid;
143//!
144//! let sg = SparseGrid::clenshaw_curtis(3, 3).unwrap();
145//! let result = sg.rule().integrate_box(
146//!     &[0.0, 0.0, 0.0], &[1.0, 1.0, 1.0],
147//!     |x| x[0] * x[1] * x[2],
148//! );
149//! assert!((result - 0.125).abs() < 1e-12);
150//! ```
151
152#[cfg(not(feature = "std"))]
153extern crate alloc;
154
155/// Implement the standard quadrature rule accessor methods.
156///
157/// Most 1-D rule types wrap a `QuadratureRule<f64>` in a field called `rule`
158/// and expose the same four accessors: `rule()`, `order()`, `nodes()`, and
159/// `weights()`. This macro generates all four with the correct attributes
160/// and documentation.
161///
162/// `$nodes_doc` customises the doc comment on `nodes()` to describe the
163/// domain (e.g. "\\[-1, 1\\]", "\\[0, ∞)", "(-∞, ∞)").
164macro_rules! impl_rule_accessors {
165    ($ty:ty, nodes_doc: $nodes_doc:expr) => {
166        impl $ty {
167            /// Returns a reference to the underlying quadrature rule.
168            #[inline]
169            #[must_use]
170            pub fn rule(&self) -> &$crate::rule::QuadratureRule<f64> {
171                &self.rule
172            }
173
174            /// Returns the number of quadrature points.
175            #[inline]
176            #[must_use]
177            pub fn order(&self) -> usize {
178                self.rule.order()
179            }
180
181            #[doc = $nodes_doc]
182            #[inline]
183            #[must_use]
184            pub fn nodes(&self) -> &[f64] {
185                &self.rule.nodes
186            }
187
188            /// Returns the weights.
189            #[inline]
190            #[must_use]
191            pub fn weights(&self) -> &[f64] {
192                &self.rule.weights
193            }
194        }
195    };
196}
197
198pub mod adaptive;
199#[cfg(feature = "std")]
200pub mod cache;
201pub mod cauchy_pv;
202pub mod clenshaw_curtis;
203pub mod cubature;
204pub mod error;
205pub mod gauss_chebyshev;
206pub mod gauss_hermite;
207pub mod gauss_jacobi;
208pub mod gauss_kronrod;
209pub mod gauss_laguerre;
210pub mod gauss_legendre;
211pub mod gauss_lobatto;
212pub mod gauss_radau;
213pub(crate) mod golub_welsch;
214pub mod oscillatory;
215pub mod result;
216pub mod rule;
217pub mod tanh_sinh;
218pub mod transforms;
219pub mod weighted;
220
221pub use adaptive::{adaptive_integrate, adaptive_integrate_with_breaks, AdaptiveIntegrator};
222pub use cauchy_pv::{pv_integrate, CauchyPV};
223pub use clenshaw_curtis::ClenshawCurtis;
224pub use error::QuadratureError;
225pub use gauss_chebyshev::{GaussChebyshevFirstKind, GaussChebyshevSecondKind};
226pub use gauss_hermite::GaussHermite;
227pub use gauss_jacobi::GaussJacobi;
228pub use gauss_kronrod::{GKPair, GaussKronrod};
229pub use gauss_laguerre::GaussLaguerre;
230pub use gauss_legendre::GaussLegendre;
231pub use gauss_lobatto::GaussLobatto;
232pub use gauss_radau::GaussRadau;
233pub use oscillatory::{
234    integrate_oscillatory_cos, integrate_oscillatory_sin, OscillatoryIntegrator, OscillatoryKernel,
235};
236pub use result::QuadratureResult;
237pub use rule::QuadratureRule;
238pub use tanh_sinh::{tanh_sinh_integrate, TanhSinh};
239pub use transforms::{
240    integrate_infinite, integrate_semi_infinite_lower, integrate_semi_infinite_upper,
241};
242pub use weighted::{weighted_integrate, WeightFunction, WeightedIntegrator};
bilby/lib.rs

bilby/
lib.rs
