//! Maths in plain Rust. Its float functions take `f16`, `f32` and `f64`, each within 4 units in the last place (ulps) of the correctly rounded result. The README lists the functions with wider bounds. The crate also has APL's tolerant comparison, exact rationals, linear algebra, number theory and polynomials. Its Python functions evaluate eagerly and return Rustymath arrays or Python scalars.
#[macro_use]
pub mod slice;
mod arc;
mod arith;
#[macro_use]
pub mod buffer;
pub mod agreement;
pub mod array;
pub mod characters;
mod complex;
pub mod convert;
mod dd;
pub mod derivative;
pub mod distribution;
pub mod domain;
pub mod error;
pub mod events;
pub mod exact;
mod exp;
pub mod facts;
mod field;
mod float;
pub mod gather;
pub mod hashing;
mod hyperbolic;
pub mod integer;
pub mod interrupt;
mod lambert;
pub mod layout;
pub mod linalg;
mod ln;
pub mod memory;
pub mod number_theory;
pub mod optimize;
pub mod parallel;
mod parts;
pub mod polynomial;
mod pow;
pub mod random;
pub mod recipe;
pub mod reduce;
pub mod search;
pub mod signal;
mod special;
mod sqrts;
pub mod stats;
pub mod structure;
mod table;
pub mod tolerant;
mod trig;
pub use arc::{RealAcos, RealAsin};
pub use arith::{
Add, BitAnd, BitOr, Bitwise, Divide, LeftShift, MathishDivide, Multiply, Nand, Negate, Nor, Not, PiTimes, Reciprocal, RightShift, Subtract, bitwise,
};
pub use complex::{Cis, divide};
pub use exp::{Exp, Expm1};
pub use float::{FUSED, Float, Native};
pub use hyperbolic::{RealAcosh, RealAtanh};
pub use ln::{Ln, Log1p, Log2, Log10, Logb, MathishLogb, RealLn};
pub use parts::{Abs, Angle, Choose, Clip, Conj, Imag, Isfinite, Isinf, Isnan, Real, Rint, Round, Scale, Sign, Trunc};
pub use pow::{Exponent, Pow, RealPow, Root};
pub use slice::{Arg, Binary, CHUNK, Flipped, KernelItem, Results, Ternary, Unary, Wrapped, pairs};
pub use special::{BetaincAt, Erfc, Erfcinv, Erfinv};
pub use sqrts::{RealSqrt, RealSqrt1mx2, RealSqrtm1mx2, RealSqrtx2m1, Sqrt, Sqrtm1mx2};
pub use trig::{Cos, Sin, Tan};
/// The crates whose types compact buffers hold, which the exported width macros name.
pub use {half, num_complex};
/// How many items a function handles between checks for an interrupt ([`interrupt::check`]), which keeps the check out of the loop over items.
pub const BETWEEN_CHECKS: usize = 1024;
/// Declares each function at the crate root, where it multiplies and adds as the target's `FUSED` says, and in `fused` and `unfused`, which choose for themselves. Each entry names the module whose function of the same name implements it, generic over `R: Native` and `const FUSED: bool`. Every wrapper is `#[inline(always)]`, because LLVM leaves a large function such as `gamma` out of line in a caller's loop, and the call stops that loop vectorizing. An entry marked `as` also declares a function object of that name, for the facts table's formula. An object with forms, such as [`Pow`] or [`Atan`], is written with `forms!` beside its function. `function_doc!(ln)` gives an entry's doc by its name. The Python package's functions take their docs from it.
macro_rules! functions {
($($(#[doc = $doc:literal])* $name:ident($($x:ident),+) in $module:ident $(as $object:ident)?;)+) => {
$(
$(#[doc = $doc])*
#[inline(always)]
pub fn $name<F: Float>($($x: F),+) -> F { functions!(@call FUSED, $module, $name, ($($x),+)) }
)+
$(functions!(@object $($object)?; $name, $module, ($($x),+));)+
/// Each function's doc by its name, as `function_doc!(ln)`, with its lines joined by newlines.
#[macro_export]
macro_rules! function_doc {
$(($name) => { concat!($($doc, "\n"),*) };)+
}
/// rustymath's functions with `mla` fused, whatever the target's `FUSED`. Code compiled with FMA enabled for itself, such as an AVX2 version of a loop chosen when the program runs, calls these on a target whose baseline lacks FMA.
pub mod fused {
use crate::Float;
$(
#[doc = concat!("[`", stringify!($name), "`](crate::", stringify!($name), "), with `mla` fused.")]
#[inline(always)]
pub fn $name<F: Float>($($x: F),+) -> F { functions!(@call true, $module, $name, ($($x),+)) }
)+
}
/// rustymath's functions with `mla` unfused, whatever the target's `FUSED`.
pub mod unfused {
use crate::Float;
$(
#[doc = concat!("[`", stringify!($name), "`](crate::", stringify!($name), "), with `mla` unfused.")]
#[inline(always)]
pub fn $name<F: Float>($($x: F),+) -> F { functions!(@call false, $module, $name, ($($x),+)) }
)+
}
};
(@object; $name:ident, $module:ident, ($($x:ident),+)) => {};
(@object $object:ident; $name:ident, $module:ident, ($($x:ident),+)) => {
#[doc = concat!("[`", stringify!($name), "`](crate::", stringify!($name), ") as a function object.")]
#[derive(Clone, Copy, Debug)]
pub struct $object;
};
(@call $fused:expr, $module:ident, $name:ident, ($($x:ident),+)) => { F::from_working(crate::$module::$name::<_, { $fused }>($($x.working()),+)) };
}
functions! {
/// e^x. `exp(∞)` is ∞, `exp(-∞)` is 0, and `exp(NaN)` is NaN.
exp(x) in exp;
/// e^x - 1, without the loss of digits near 0 that `exp(x) - 1` has.
expm1(x) in exp;
/// The natural logarithm of x. `ln(0)` is -∞, `ln(∞)` is ∞, and a negative real x or NaN gives NaN.
ln(x) in ln;
/// ln(1 + x), without the loss of digits near 0 that `ln(1 + x)` has. `log1p(-1)` is -∞.
log1p(x) in ln;
/// The logarithm of x to base 2. It is k exactly at 2^k.
log2(x) in ln;
/// The logarithm of x to base 10. It is k exactly at 10^k, and at the float nearest 10^k where 10^k isn't a float.
log10(x) in ln;
/// The logarithm of x to base b, ln x / ln b. Base 2 and base 10 take `log2` and `log10`, which are exact at powers of their base.
logb(b, x) in ln;
/// sin x. `sin(±0)` is ±0, and `sin(±∞)` and `sin(NaN)` are NaN. Beyond a magnitude of 2^23 at f64 or 2^20 at f32, x takes a slower, exact reduction.
sin(x) in trig;
/// cos x. `cos(±0)` is 1, and `cos(±∞)` and `cos(NaN)` are NaN. Beyond the magnitude that `sin` gives, x takes the same exact reduction.
cos(x) in trig;
/// tan x. `tan(±0)` is ±0, and `tan(±∞)` and `tan(NaN)` are NaN. Beyond the magnitude that `sin` gives, x takes the same exact reduction.
tan(x) in trig;
/// asin x. A real x gives a result from -π/2 to π/2 and NaN beyond a magnitude of 1.
asin(x) in arc as Asin;
/// acos x. A real x gives a result from 0 to π and NaN beyond a magnitude of 1.
acos(x) in arc as Acos;
/// atan x. A real x gives a result from -π/2 to π/2.
atan(x) in arc as Atan;
/// The angle of the point (x, y) from the positive x axis, from -π to π, with y first as in C's `atan2`. Zeros and infinities give C99's angles: atan2(±0, -0) is ±π, and atan2(±∞, ±∞) is an odd multiple of ±π/4.
atan2(y, x) in arc as Atan2;
/// sinh x. It stays finite up to where sinh itself overflows, beyond the point where e^x does.
sinh(x) in hyperbolic as Sinh;
/// cosh x. It stays finite up to where cosh itself overflows, beyond the point where e^x does.
cosh(x) in hyperbolic as Cosh;
/// tanh x.
tanh(x) in hyperbolic as Tanh;
/// asinh x.
asinh(x) in hyperbolic as Asinh;
/// acosh x. A real x below 1 gives NaN.
acosh(x) in hyperbolic as Acosh;
/// atanh x. ±1 gives ±∞, and a real x of magnitude beyond 1 gives NaN.
atanh(x) in hyperbolic as Atanh;
/// x^y. Special values follow IEEE 754's `pow`: x^0 and 1^y are 1, even for NaN, and a negative x has a power only for a whole y.
pow(x, y) in pow;
/// The real nth root of y, computed without rounding 1/n first. A negative y has one only for an odd whole n. n = 0, an infinite n and NaN give NaN.
root(n, y) in pow;
/// The cube root of x.
cbrt(x) in pow as Cbrt;
/// √(x² + y²), with no overflow or underflow in between. An infinite argument gives ∞, even beside NaN.
hypot(x, y) in pow as Hypot;
/// √(1 - x²), without the digits that squaring x first loses near |x| = 1. A real x of magnitude beyond 1 gives NaN.
sqrt1mx2(x) in sqrts as Sqrt1mx2;
/// √(1 + x²), with no overflow where x² would overflow.
sqrt1px2(x) in sqrts as Sqrt1px2;
/// √(x - 1) √(x + 1): √(x² - 1) with the sign of a real x, without the digits that squaring x first loses near |x| = 1. A real x of magnitude below 1 gives NaN.
sqrtx2m1(x) in sqrts as Sqrtx2m1;
/// Γ(x). Γ(±0) is ±∞, a negative whole x and -∞ give NaN, and Γ(∞) is ∞.
gamma(x) in special as Gamma;
/// Γ(x + 1), which is x! for a whole x from 0.
factorial(x) in special as Factorial;
/// ln |Γ(x)|. It is ∞ at 0, at the negative whole numbers and at ±∞. Near its zeros at 1 and 2 its error is within 4 ulps of 2, not of the result.
lgamma(x) in special as Lgamma;
/// ψ(x), the digamma function: the derivative of ln Γ(x), Γ′(x)/Γ(x). ψ(+0) is -∞, ψ(-0) is ∞, a negative whole x and -∞ give NaN, and ψ(∞) is ∞. Near its zeros its error is within 4 ulps of 1, not of the result, and below -20 within 4 ulps of ln(1 - x).
digamma(x) in special as Digamma;
/// ψ⁽ⁿ⁾(x), the polygamma function: the nth derivative of ψ, for a whole n from 0 to 170. For n from 1, x from 0 down gives NaN, as do a fractional or negative n and an n above 170.
polygamma(n, x) in special as Polygamma;
/// B(a, b) = Γ(a) Γ(b)/Γ(a + b), the beta function. A pole of Γ(a) or Γ(b) gives NaN, and a pole of Γ(a + b) alone gives 0.
beta(a, b) in special as Beta;
/// C(n, k), the binomial coefficient. Whole n and k from 0 multiply out exactly while the result times |n| + k is a whole number the type holds, and a negative whole n gives (-1)^k C(k - n - 1, k). Other n with a whole k from 0 to 64 multiply out n(n - 1)…(n - k + 1)/k!, within 2k ulps. Other real arguments take Γ(n + 1)/(Γ(k + 1) Γ(n - k + 1)) in double words, within 4 ulps in tests against mpmath, for large n and negative arguments too.
binom(n, k) in special as Binom;
/// ∂C(n, k)/∂n, the binomial coefficient's slope in n. For a whole k from 0 to 64 it is the slope of the polynomial n(n - 1)…(n - k + 1)/k! by the product rule, within 2k ulps. Other arguments take C(n, k)(ψ(n + 1) - ψ(n - k + 1)), with the ψ difference in double words where both arguments are from 10 up, with the polynomial's slope where a whole n below a whole k puts those ψ at a pole. Other poles of ψ give NaN.
binom_dn(n, k) in special as BinomDn;
/// ∂C(n, k)/∂k, the binomial coefficient's slope in k: C(n, k)(ψ(n - k + 1) - ψ(k + 1)). At a whole n from 0 below a whole k, it is the negative of `binom_dn`. Other poles of ψ, as at a negative whole k, give NaN. For a whole k from 0 to 64, its error is within 2k ulps of C(n, k) times ψ(n - k + 1) - ψ(k + 1), plus the error of that difference, so near n = 2k, where the two ψ cancel, its relative error has no bound.
binom_dk(n, k) in special as BinomDk;
/// The error function, 2/√π times the integral of e^-t² from 0 to x.
erf(x) in special as Erf;
/// 1 - erf(x), without the loss of digits for large x that `1 - erf(x)` has.
erfc(x) in special;
/// The x with erf(x) = p. ±1 gives ±∞, and a magnitude beyond 1 gives NaN.
erfinv(p) in special;
/// The x with erfc(x) = q. 0 gives ∞, 2 gives -∞, and q outside 0 to 2 gives NaN.
erfcinv(q) in special;
/// W0(x), the principal branch of the Lambert W function: the w of at least -1 with w e^w = x. A real x below -1/e gives NaN, except within half an ulp of -1/e, where W0 is -1.
lambertw(x) in lambert as LambertW;
/// The regularized lower incomplete gamma function P(a, x): the integral of t^(a-1) e^-t from 0 to x, divided by Γ(a). a must be positive and x nonnegative, or the result is NaN.
gammainc(a, x) in special as Gammainc;
/// The regularized upper incomplete gamma function Q(a, x) = 1 - P(a, x), without the loss of digits that `1 - gammainc(a, x)` has where Q is small.
gammaincc(a, x) in special as Gammaincc;
/// The x with P(a, x) = p, for p from 0 to 1.
gammaincinv(a, p) in special as Gammaincinv;
/// The regularized incomplete beta function I_x(a, b): the integral of t^(a-1) (1-t)^(b-1) from 0 to x, divided by B(a, b). a and b must be positive and x from 0 to 1, or the result is NaN.
betainc(a, b, x) in special as Betainc;
/// The x with I_x(a, b) = p, for p from 0 to 1.
betaincinv(a, b, p) in special as Betaincinv;
}