r2rs-nmath 0.1.1

//! # Special Functions of Mathematics
//!
//! Special mathematical functions related to the beta and gamma functions.
//!
//! ## Arguments
//!
//! * a, b: non-negative numeric
//! * x, n: numeric
//! * k, deriv: integer
//!
//! ## Details
//!
//! The functions beta and lbeta return the beta function and the natural logarithm of the beta
//! function,
//!
//! $B(a,b) = \Gamma(a)\Gamma(b)/\Gamma(a+b)$.
//!
//! The formal definition is
//!
//! $\int_0^1 t^{a-1} (1-t)^{b-1} dt$
//!
//! (Abramowitz and Stegun section 6.2.1, page 258). Note that it is only defined in R for
//! non-negative a and b, and is infinite if either is zero.
//!
//! The functions gamma and lgamma return the gamma function Γ(x) and the natural logarithm of
//! the absolute value of the gamma function. The gamma function is defined by (Abramowitz and Stegun section 6.1.1, page 255)
//!
//! $\gamma(x) = \int_0^\infty t^{x-1} \text{exp}(-t) dt$
//!
//! for all real x except zero and negative integers (when NaN is returned). There will be a
//! warning on possible loss of precision for values which are too close (within about 1e-8)
//! to a negative integer less than -10.
//!
//! factorial(x) (x! for non-negative integer x) is defined to be gamma(x+1) and lfactorial to be
//! lgamma(x+1).
//!
//! The functions digamma and trigamma return the first and second derivatives of the logarithm
//! of the gamma function. psigamma(x, deriv) (deriv >= 0) computes the deriv-th derivative of $\psi(x)$.
//!
//! $digamma(x) = \psi(x) = \frac{d}{dx}\{ln \Gamma(x)\} = \Gamma'(x) / \Gamma(x)$
//!
//! \psi and its derivatives, the psigamma() functions, are often called the ‘polygamma’ functions,
//! e.g. in Abramowitz and Stegun (section 6.4.1, page 260); and higher derivatives (deriv = 2:4)
//! have occasionally been called ‘tetragamma’, ‘pentagamma’, and ‘hexagamma’.
//!
//! The functions choose and lchoose return binomial coefficients and the logarithms of their
//! absolute values. Note that choose(n, k) is defined for all real numbers n and integer k.
//! For $k ≥ 1$ it is defined as $n(n-1)…(n-k+1) / k!$, as 1 for $k = 0$ and as 0 for negative k.
//! Non-integer values of k are rounded to an integer, with a warning.
//!
//! choose(*, k) uses direct arithmetic (instead of \[l\]gamma calls) for small k, for speed and
//! accuracy reasons. Note the function combn (package utils) for enumeration of all possible
//! combinations.
//!
//! The gamma, lgamma, digamma and trigamma functions are internal generic primitive functions:
//! methods can be defined for them individually or via the Math group generic.
//!
//! ## Source
//!
//! gamma, lgamma, beta and lbeta are based on C translations of Fortran subroutines by W.
//! Fullerton of Los Alamos Scientific Laboratory (now available as part of SLATEC).
//!
//! digamma, trigamma and psigamma are based on
//!
//! Amos, D. E. (1983). A portable Fortran subroutine for derivatives of the psi function,
//! Algorithm 610, ACM Transactions on Mathematical Software 9(4), 494–502.
//!
//! ## References
//!
//! Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) The New S Language. Wadsworth &
//! Brooks/Cole. (For gamma and lgamma.)
//!
//! Abramowitz, M. and Stegun, I. A. (1972) Handbook of Mathematical Functions. New York: Dover.
//! https://!en.wikipedia.org/wiki/Abramowitz_and_Stegun provides links to the full text which is
//! in public domain.
//! Chapter 6: Gamma and Related Functions.
//!
//! ## See Also
//!
//! Arithmetic for simple, sqrt for miscellaneous mathematical functions and Bessel for the real
//! Bessel functions.
//!
//! For the incomplete gamma function see pgamma.
//!
//! ## Examples
//!
//! Basic Usage:
//! ```rust
//! # use r2rs_nmath::func::*;
//! # use strafe_type::FloatConstraint;
//! # use std::{fs::File, io::Write};
//! println!("{}", choose(5, 2));
//! println!("{:e}", factorial(100).unwrap());
//! println!("{}", lfactorial(10000));
//! # let mut f = File::create("src/func/betagam/doctest_out/betagam.md").unwrap();
//! # writeln!(f, "```output").unwrap();
//! # writeln!(f, "{}", choose(5, 2)).unwrap();
//! # writeln!(f, "{:e}", factorial(100).unwrap()).unwrap();
//! # writeln!(f, "{}", lfactorial(10000)).unwrap();
//! # writeln!(f, "```").unwrap();
//! ```
#![cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = include_str!("doctest_out/betagam.md")))]
//!
//! Pascal's Triangle:
//! ```rust
//! # use r2rs_nmath::func::choose;
//! # use strafe_type::FloatConstraint;
//! # use std::{fs::File, io::Write};
//! for n in 0..=10 {
//!     for i in 0..=n {
//!         print!("{} ", choose(n, i))
//!     }
//!     println!();
//! }
//! # let mut f = File::create("src/func/betagam/doctest_out/pascal.md").unwrap();
//! # writeln!(f, "```output").unwrap();
//! # for n in 0..=10 {
//! #    for i in 0..=n {
//! #        write!(f, "{} ", choose(n, i)).unwrap();
//! #    }
//! #    writeln!(f, "").unwrap();
//! # }
//! # writeln!(f, "```").unwrap();
//! ```
#![cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = include_str!("doctest_out/pascal.md")))]
//!
//! Gamma Plot:
//! ```rust
//! // Gamma plot
//! # use r2rs_base::traits::StatisticalSlice;
//! # use r2rs_nmath::func::gamma;
//! # use strafe_plot::prelude::*;
//! # use strafe_type::FloatConstraint;
//! # use std::fs::rename;
//! let mut plot = Plot::new();
//! plot.with_options(PlotOptions {
//!     y_min: Some(-20.0),
//!     y_max: Some(20.0),
//!     title: "Γ(x)".to_string(),
//!     x_axis_label: "x".to_string(),
//!     y_axis_label: "gamma(x)".to_string(),
//!     ..Default::default()
//! });
//!
//! // gamma has 1st order poles at 0, -1, -2, ...
//! for sections in [-3.0, -2.0, -1.0, 0.0, 4.0].windows(2) {
//!     let left = sections[0];
//!     let right = sections[1];
//!
//!     let x = <[f64]>::sequence(left + 1e-4, right - 1e-4, 200);
//!     let y = x
//!         .iter()
//!         .map(|x| gamma(x).unwrap().unwrap())
//!         .collect::<Vec<_>>();
//!
//!     plot.with_plottable(Line {
//!         x,
//!         y,
//!         color: RED,
//!         ..Default::default()
//!     });
//! }
//!
//! plot.with_plottable(HorizontalLine {
//!     y: 0.0,
//!     color: BLUE,
//!     size: 1,
//!     dash: true,
//!     ..Default::default()
//! })
//! .with_plottable(VerticalLine {
//!     x: -2.0,
//!     color: BLUE,
//!     size: 1,
//!     dash: true,
//!     ..Default::default()
//! })
//! .with_plottable(VerticalLine {
//!     x: -1.0,
//!     color: BLUE,
//!     size: 1,
//!     dash: true,
//!     ..Default::default()
//! })
//! .with_plottable(VerticalLine {
//!     x: 0.0,
//!     color: BLUE,
//!     size: 1,
//!     dash: true,
//!     ..Default::default()
//! })
//! .plot(&SVGBackend::new("gamma.svg", (1024, 768)).into_drawing_area())
//! .unwrap();
//! # rename(
//! #     format!("gamma.svg"),
//! #     format!("src/func/betagam/doctest_out/gamma.svg"),
//! # )
//! # .unwrap();
//! ```
#![cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = embed_doc_image::embed_image!("gamma", "src/func/betagam/doctest_out/gamma.svg")))]
#![cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = "![Gamma Plot][gamma]"))]
//!
//! Psigamma Plot:
//! ```rust
//! # use r2rs_base::traits::StatisticalSlice;
//! # use r2rs_nmath::func::{gamma, lgamma, psigamma};
//! # use strafe_plot::prelude::{*, full_palette::GREY_600};
//! # use strafe_type::FloatConstraint;
//! # use std::fs::rename;
//! let x = <[f64]>::sequence(0.1, 4.0, 201);
//!
//! let dx = x.windows(2).map(|xs| xs[1] - xs[0]).collect::<Vec<_>>();
//! let x_temp = x[0..201].to_vec();
//!
//! let y = x
//!     .iter()
//!     .map(|x| gamma(x).unwrap().unwrap())
//!     .collect::<Vec<_>>();
//!
//! let root = SVGBackend::new("psigamma.svg", (1024, 768)).into_drawing_area();
//!
//! Plot::new()
//!     .with_options(PlotOptions {
//!         plot_right: 1.0 / 3.0,
//!         plot_bottom: 0.5,
//!         title: "gamma".to_string(),
//!         x_axis_label: "x".to_string(),
//!         y_axis_label: "y".to_string(),
//!         ..Default::default()
//!     })
//!     .with_plottable(HorizontalLine {
//!         y: 0.0,
//!         color: GREY_600,
//!         size: 1,
//!         ..Default::default()
//!     })
//!     .with_plottable(Line {
//!         x: x.clone(),
//!         y,
//!         color: RED,
//!         ..Default::default()
//!     })
//!     .plot(&root)
//!     .unwrap();
//!
//! let y = x
//!     .iter()
//!     .map(|x| lgamma(x).unwrap().unwrap())
//!     .collect::<Vec<_>>();
//!
//! Plot::new()
//!     .with_options(PlotOptions {
//!         plot_left: 1.0 / 3.0,
//!         plot_right: 2.0 / 3.0,
//!         plot_bottom: 0.5,
//!         title: "lgamma".to_string(),
//!         x_axis_label: "x".to_string(),
//!         y_axis_label: "y".to_string(),
//!         ..Default::default()
//!     })
//!     .with_plottable(HorizontalLine {
//!         y: 0.0,
//!         color: GREY_600,
//!         size: 1,
//!         ..Default::default()
//!     })
//!     .with_plottable(Line {
//!         x: x.clone(),
//!         y: y.clone(),
//!         color: RED,
//!         ..Default::default()
//!     })
//!     .plot(&root)
//!     .unwrap();
//!
//! let dy = y
//!     .windows(2)
//!     .zip(dx.iter())
//!     .map(|(ys, dx)| (ys[1] - ys[0]) / dx)
//!     .collect::<Vec<_>>();
//! let y = x
//!     .iter()
//!     .map(|x| psigamma(x, 0).unwrap().unwrap())
//!     .collect::<Vec<_>>();
//!
//! Plot::new()
//!     .with_options(PlotOptions {
//!         plot_left: 2.0 / 3.0,
//!         plot_right: 1.0,
//!         plot_bottom: 0.5,
//!         title: "digamma".to_string(),
//!         x_axis_label: "x".to_string(),
//!         y_axis_label: "y".to_string(),
//!         ..Default::default()
//!     })
//!     .with_plottable(HorizontalLine {
//!         y: 0.0,
//!         color: GREY_600,
//!         size: 1,
//!         ..Default::default()
//!     })
//!     .with_plottable(Line {
//!         x: x.clone(),
//!         y: y.clone(),
//!         color: RED,
//!         ..Default::default()
//!     })
//!     .with_plottable(Line {
//!         x: x_temp.clone(),
//!         y: dy,
//!         color: BLUE,
//!         dash: true,
//!         ..Default::default()
//!     })
//!     .plot(&root)
//!     .unwrap();
//!
//! let dy = y
//!     .windows(2)
//!     .zip(dx.iter())
//!     .map(|(ys, dx)| (ys[1] - ys[0]) / dx)
//!     .collect::<Vec<_>>();
//! let y = x
//!     .iter()
//!     .map(|x| psigamma(x, 1).unwrap().unwrap())
//!     .collect::<Vec<_>>();
//!
//! Plot::new()
//!     .with_options(PlotOptions {
//!         plot_right: 1.0 / 3.0,
//!         plot_top: 0.5,
//!         plot_bottom: 1.0,
//!         title: "trigamma".to_string(),
//!         x_axis_label: "x".to_string(),
//!         y_axis_label: "y".to_string(),
//!         ..Default::default()
//!     })
//!     .with_plottable(HorizontalLine {
//!         y: 0.0,
//!         color: GREY_600,
//!         size: 1,
//!         ..Default::default()
//!     })
//!     .with_plottable(Line {
//!         x: x.clone(),
//!         y: y.clone(),
//!         color: RED,
//!         ..Default::default()
//!     })
//!     .with_plottable(Line {
//!         x: x_temp.clone(),
//!         y: dy,
//!         color: BLUE,
//!         dash: true,
//!         ..Default::default()
//!     })
//!     .plot(&root)
//!     .unwrap();
//!
//! let dy = y
//!     .windows(2)
//!     .zip(dx.iter())
//!     .map(|(ys, dx)| (ys[1] - ys[0]) / dx)
//!     .collect::<Vec<_>>();
//! let y = x
//!     .iter()
//!     .map(|x| psigamma(x, 2).unwrap().unwrap())
//!     .collect::<Vec<_>>();
//!
//! Plot::new()
//!     .with_options(PlotOptions {
//!         plot_left: 1.0 / 3.0,
//!         plot_right: 2.0 / 3.0,
//!         plot_top: 0.5,
//!         plot_bottom: 1.0,
//!         title: "psigamma(deriv=2)".to_string(),
//!         x_axis_label: "x".to_string(),
//!         y_axis_label: "y".to_string(),
//!         ..Default::default()
//!     })
//!     .with_plottable(HorizontalLine {
//!         y: 0.0,
//!         color: GREY_600,
//!         size: 1,
//!         ..Default::default()
//!     })
//!     .with_plottable(Line {
//!         x: x.clone(),
//!         y: y.clone(),
//!         color: RED,
//!         ..Default::default()
//!     })
//!     .with_plottable(Line {
//!         x: x_temp.clone(),
//!         y: dy,
//!         color: BLUE,
//!         dash: true,
//!         ..Default::default()
//!     })
//!     .plot(&root)
//!     .unwrap();
//!
//! let dy = y
//!     .windows(2)
//!     .zip(dx.iter())
//!     .map(|(ys, dx)| (ys[1] - ys[0]) / dx)
//!     .collect::<Vec<_>>();
//! let y = x
//!     .iter()
//!     .map(|x| psigamma(x, 3).unwrap().unwrap())
//!     .collect::<Vec<_>>();
//!
//! Plot::new()
//!     .with_options(PlotOptions {
//!         plot_left: 2.0 / 3.0,
//!         plot_right: 1.0,
//!         plot_top: 0.5,
//!         plot_bottom: 1.0,
//!         title: "psigamma(deriv=3)".to_string(),
//!         x_axis_label: "x".to_string(),
//!         y_axis_label: "y".to_string(),
//!         ..Default::default()
//!     })
//!     .with_plottable(HorizontalLine {
//!         y: 0.0,
//!         color: GREY_600,
//!         size: 1,
//!         ..Default::default()
//!     })
//!     .with_plottable(Line {
//!         x: x.clone(),
//!         y,
//!         color: RED,
//!         ..Default::default()
//!     })
//!     .with_plottable(Line {
//!         x: x_temp.clone(),
//!         y: dy,
//!         color: BLUE,
//!         dash: true,
//!         ..Default::default()
//!     })
//!     .plot(&root)
//!     .unwrap();
//! # drop(root);
//! # rename(
//! #     format!("psigamma.svg"),
//! #     format!("src/func/betagam/doctest_out/psigamma.svg"),
//! # )
//! # .unwrap();
//! ```
#![cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = embed_doc_image::embed_image!("psigamma", "src/func/betagam/doctest_out/psigamma.svg")))]
#![cfg_attr(
    feature = "doc_outputs",
    cfg_attr(all(), doc = "![Psigamma Plot][psigamma]")
)]
//!
//! Binomial theorem for $n = \frac{1}{2}$
//!
//! $\sqrt{1+x} = (1+x)^{\frac{1}{2}} = \sum_{k=0}^\infty \text{choose}(\frac{1}{2}, k) * x^k$
//! ```rust
//! # use r2rs_nmath::func::choose;
//! # use strafe_type::FloatConstraint;
//! let k = 0..10; // 10 is sufficient for ~ 9 digit precision:
//! println!("{}", 1.25_f64.sqrt());
//! let sum = k
//!     .map(|k| choose(1.0 / 2.0, k).unwrap() * 0.25_f64.powi(k))
//!     .sum::<f64>();
//! println!("{sum}");
//! # use std::{fs::File, io::Write};
//! # let mut f = File::create("src/func/betagam/doctest_out/binomial_theorem.md").unwrap();
//! # writeln!(f, "```output").unwrap();
//! # writeln!(f, "{}", 1.25_f64.sqrt()).unwrap();
//! # writeln!(f, "{sum}").unwrap();
//! # writeln!(f, "```").unwrap();
//! ```
#![cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = include_str!("doctest_out/binomial_theorem.md")))]

mod beta_;
mod choose_;
mod gamma_;
mod gamma_cody_;
mod lbeta_;
mod lgamma_;
mod lgammacor_;
mod polygamma_;

pub use self::{
    beta_::*, choose_::*, factorial::*, gamma_::*, gamma_cody_::*, lbeta_::*, lgamma_::*,
    lgammacor_::*, polygamma_::*,
};

mod factorial;

#[cfg(test)]
mod tests;

#[cfg(all(test, feature = "enable_proptest"))]
mod proptests;

#[cfg(all(test, feature = "enable_covtest"))]
mod covtests;