gauss-quad 0.3.1

// Copyright 2019-2024 Dominique Dresen
// Copyright 2023-2026 Johanna Sörngård
// SPDX-License-Identifier: MIT OR Apache-2.0

//! # gauss-quad
//!
//! **gauss-quad** is a [Gaussian quadrature](https://en.wikipedia.org/wiki/Gaussian_quadrature) library for numerical integration.
//!
//! ## Quadrature rules
//!
//! **gauss-quad** implements the following quadrature rules:
//! * [Gauss-Legendre](https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_quadrature)
//! * [Gauss-Jacobi](https://en.wikipedia.org/wiki/Gauss%E2%80%93Jacobi_quadrature)
//! * [Gauss-Laguerre](https://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature) (generalized)
//! * [Gauss-Hermite](https://en.wikipedia.org/wiki/Gauss%E2%80%93Hermite_quadrature)
//! * [Gauss-Chebyshev](https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature)
//! * [Trapezoid](https://en.wikipedia.org/wiki/Trapezoidal_rule)
//! * [Midpoint](https://en.wikipedia.org/wiki/Riemann_sum#Midpoint_rule)
//! * [Simpson](https://en.wikipedia.org/wiki/Simpson%27s_rule)
//!
//! ## Using **gauss-quad**
//!
//! To use any of the quadrature rules in your project, first initialize the rule with
//! a specified degree and then you can use it for integration. The degree
//! is a non-zero positive integer that determines the number of nodes used in the quadrature rule,
//! and as such the number of points at which the integrand is evaluated.
//!
//! ```
//! use gauss_quad::GaussLegendre;
//! // This macro is used in these docs to compare floats.
//! // The assertion succeeds if the two sides are within floating point error,
//! // or an optional epsilon.
//! use approx::assert_abs_diff_eq;
//! use core::num::NonZeroUsize;
//!
//! // initialize the quadrature rule
//! let degree = NonZeroUsize::new(10).unwrap();
//! let quad = GaussLegendre::new(degree);
//!
//! // Use the rule to integrate a function
//! let left_bound = 0.0;
//! let right_bound = 1.0;
//! let integral = quad.integrate(left_bound, right_bound, |x| x * x);
//! assert_abs_diff_eq!(integral, 1.0 / 3.0);
//! ```
//!
//! Select the degree, n, such that 2n-1 is the largest degree of polynomial that you want to integrate with the rule.
//!
//! ## Setting up a quadrature rule
//!
//! Using a quadrature rule takes two steps:
//! 1. Initialization
//! 2. Integration
//!
//! First, rules must be initialized using some specific input parameters.
//!
//! Then, you can integrate functions using those rules:
//!
//! ```
//! # use gauss_quad::*;
//! # let degree = core::num::NonZeroUsize::new(5).unwrap();
//! # let alpha = FiniteAboveNegOneF64::new(1.2).unwrap();
//! # let beta = FiniteAboveNegOneF64::new(1.2).unwrap();
//! # let a = 0.0;
//! # let b = 1.0;
//! # let c = -10.;
//! # let d = 100.;
//! let gauss_legendre = GaussLegendre::new(degree);
//! // Integrate on the domain [a, b]
//! let x_cubed = gauss_legendre.integrate(a, b, |x| x * x * x);
//!
//! let gauss_jacobi = GaussJacobi::new(degree, alpha, beta);
//! // Integrate on the domain [c, d]
//! let double_x = gauss_jacobi.integrate(c, d, |x| 2.0 * x);
//!
//! let gauss_laguerre = GaussLaguerre::new(degree, alpha);
//! // No explicit domain, Gauss-Laguerre integration is done on the domain [0, ∞).
//! let piecewise = gauss_laguerre.integrate(|x| if x > 0.0 && x < 2.0 { x } else { 0.0 });
//!
//! let gauss_hermite = GaussHermite::new(degree);
//! // Again, no explicit domain since Gauss-Hermite integration is done over the domain (-∞, ∞).
//! let golden_polynomial = gauss_hermite.integrate(|x| x * x - x - 1.0);
//! ```
//!
//! ## Different rules have different requirements
//!
//! Quadrature rules are only defined for a certain set of input values.
//! They take parameters of types that guarantee valid input.
//!
//! [`GaussJacobi`] for example, requires two parameters named "alpha" and "beta",
//! and they must be larger than -1.0. Therefore they are of type [`FiniteAboveNegOneF64`],
//! which can only be created from finite non-NAN values above -1.0:
//!
//! ```
//! use gauss_quad::{GaussJacobi, FiniteAboveNegOneF64};
//! use core::num::NonZeroUsize;
//!
//! let degree = NonZeroUsize::new(10).unwrap();
//! let alpha = FiniteAboveNegOneF64::new(0.1).unwrap();
//!
//! // Trying to create a `beta` below -1.0 results in `None`.
//! let beta = FiniteAboveNegOneF64::new(-1.1);
//! assert!(beta.is_none());
//!
//! // This is valid:
//! let beta = FiniteAboveNegOneF64::new(-0.5).unwrap();
//! let quad = GaussJacobi::new(degree, alpha, beta);
//! ```
//!
//! [`GaussLaguerre`] is used to evaluate an improper integral over the domain [0, ∞).
//! This means no domain bounds are needed in the `integrate` call.
//!
//! ```
//! # use gauss_quad::{GaussLaguerre, FiniteAboveNegOneF64};
//! # use approx::assert_abs_diff_eq;
//! # use core::f64::consts::PI;
//! # use core::num::NonZeroUsize;
//! // Initialize the quadrature rule
//! let degree = NonZeroUsize::new(2).unwrap();
//! let alpha = FiniteAboveNegOneF64::new(0.5).unwrap();
//! let quad = GaussLaguerre::new(degree, alpha);
//!
//! // Use the rule to integrate a function, note the lack of domain bounds
//! let integral = quad.integrate(|x| x * x);
//!
//! assert_abs_diff_eq!(integral, 15.0 * PI.sqrt() / 8.0, epsilon = 1e-14);
//! ```
//!
//! Make sure to read the specific quadrature rule's documentation before using it.
//!
//! ## Passing functions to quadrature rules
//!
//! The `integrate` method takes any integrand that implements the [`FnMut(f64) -> f64`](FnMut) trait, i.e. functions of
//! one `f64` parameter that can modify internal state.
//!
//! ```
//! # use gauss_quad::legendre::GaussLegendre;
//! # use approx::assert_abs_diff_eq;
//! # use core::num::NonZeroUsize;
//!
//! // Initialize the quadrature rule
//! let degree = NonZeroUsize::new(2).unwrap();
//! let quad = GaussLegendre::new(degree);
//!
//! // Use the rule to integrate a closure
//! let left_bound = 0.0;
//! let right_bound = 1.0;
//!
//! let integral = quad.integrate(left_bound, right_bound, |x| x * x);
//!
//! assert_abs_diff_eq!(integral, 1.0 / 3.0);
//!
//! // You can also pass a function pointer
//! fn x_cubed(x: f64) -> f64 {
//!    x * x * x
//! }
//!
//! let integral_x_cubed = quad.integrate(left_bound, right_bound, x_cubed);
//! assert_abs_diff_eq!(integral_x_cubed, 1.0 / 4.0);
//! ```
//! The parallel `par_integrate` methods instead take a function that can't modify any state, [`Fn(f64) -> f64`](Fn), and is [`Sync`] to enable parallel evaluation.
//!
//! ## Double integrals
//!
//! It is possible to use this crate to do double and higher integrals:
//!
//! ```
//! # use gauss_quad::legendre::GaussLegendre;
//! # use approx::assert_abs_diff_eq;
//! let rule = GaussLegendre::new(3.try_into().unwrap());
//!
//! // integrate x^2*y over the triangle in the xy-plane where x ϵ [0, 1] and y ϵ [0, x]:
//! let double_int = rule.integrate(0.0, 1.0, |x| rule.integrate(0.0, x, |y| x * x * y));
//!
//! assert_abs_diff_eq!(double_int, 0.1);
//! ```
//!
//! However, the time complexity of the integration then scales with the number of nodes to
//! the power of the depth of the integral, e.g. O(n³) for triple integrals.
//!
//! ## Feature flags
//!
//! `serde`: implements the [`Serialize`](serde::Serialize) and [`Deserialize`](serde::Deserialize) traits from
//! the [`serde`] crate for the quadrature rule structs.
//!
//! `rkyv`: implements the [`Serialize`](rkyv::Serialize), [`Deserialize`](rkyv::Deserialize), and [`Archive`](rkyv::Archive) traits from the [`rkyv`] crate.
//!
//! `rayon`: enables a parallel version of the `integrate` function on the quadrature rule structs. Can speed up integration if evaluating the integrand is expensive (takes ≫100 µs).
//! The [`rayon`] crate depends on the standard library, so this also enables the `std` feature.
//!
//! `zerocopy`: imlements the [`KnownLayout`](zerocopy::KnownLayout) trait from the [`zerocopy`] crate for the quadrature rule structs.
//!
//! One of the below features must be enabled:
//!
//! `libm` (*enabled by default*): depends on the [`libm`](https://docs.rs/libm/latest/libm/) crate and uses it as the math backend.
//! Does nothing if the `std` feature is enabled.
//!
//! `std`: links the standard library and uses it as the math backend.
//! When this feature is disabled the crate is `no_std` compatible.

// Only enable the nightly `doc_cfg` feature when
// the `docsrs` configuration attribute is defined.
// The config in Cargo.toml means that this is defined when we are on docs.rs (which uses the nightly compiler)
// or if the environment variable RUSTDOCFLAGS is set as `RUSTDOCFLAGS="--cfg docsrs"`.
// This lets us get a tag on docs.rs that says "Available on crate feature ... only".
#![cfg_attr(docsrs, feature(doc_cfg))]
#![forbid(unsafe_code)]
#![forbid(clippy::unwrap_used)]
#![forbid(clippy::expect_used)]
#![no_std]

extern crate alloc;
#[cfg(feature = "std")]
extern crate std;

pub mod chebyshev;
mod data_api;
mod golub_welsch;
pub mod hermite;
pub mod jacobi;
pub mod laguerre;
pub mod legendre;
mod math;
pub mod midpoint;
pub mod simpson;
pub mod trapezoid;

#[doc(inline)]
pub use chebyshev::{GaussChebyshevFirstKind, GaussChebyshevSecondKind};
#[doc(inline)]
pub use data_api::{
    FiniteAboveNegOneF64, InfNanNegOneOrLessError, Node, ParseFiniteAboveNegOneF64Error, Weight,
};
#[doc(inline)]
pub use hermite::GaussHermite;
#[doc(inline)]
pub use jacobi::GaussJacobi;
#[doc(inline)]
pub use laguerre::GaussLaguerre;
#[doc(inline)]
pub use legendre::GaussLegendre;
#[doc(inline)]
pub use midpoint::Midpoint;
#[doc(inline)]
pub use simpson::Simpson;
#[doc(inline)]
pub use trapezoid::Trapezoid;