rustymath_core/lib.rs
1//! 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.
2#[macro_use]
3pub mod slice;
4mod arc;
5mod arith;
6#[macro_use]
7pub mod buffer;
8pub mod agreement;
9pub mod array;
10pub mod characters;
11mod complex;
12pub mod convert;
13mod dd;
14pub mod derivative;
15pub mod distribution;
16pub mod domain;
17pub mod error;
18pub mod events;
19pub mod exact;
20mod exp;
21pub mod facts;
22mod field;
23mod float;
24pub mod gather;
25pub mod hashing;
26mod hyperbolic;
27pub mod integer;
28pub mod interrupt;
29mod lambert;
30pub mod layout;
31pub mod linalg;
32mod ln;
33pub mod memory;
34pub mod number_theory;
35pub mod optimize;
36pub mod parallel;
37mod parts;
38pub mod polynomial;
39mod pow;
40pub mod random;
41pub mod recipe;
42pub mod reduce;
43pub mod search;
44pub mod signal;
45mod special;
46mod sqrts;
47pub mod stats;
48pub mod structure;
49mod table;
50pub mod tolerant;
51mod trig;
52
53pub use arc::{RealAcos, RealAsin};
54pub use arith::{
55 Add, BitAnd, BitOr, Bitwise, Divide, LeftShift, MathishDivide, Multiply, Nand, Negate, Nor, Not, PiTimes, Reciprocal, RightShift, Subtract, bitwise,
56};
57pub use complex::{Cis, divide};
58pub use exp::{Exp, Expm1};
59pub use float::{FUSED, Float, Native};
60pub use hyperbolic::{RealAcosh, RealAtanh};
61pub use ln::{Ln, Log1p, Log2, Log10, Logb, MathishLogb, RealLn};
62pub use parts::{Abs, Angle, Choose, Clip, Conj, Imag, Isfinite, Isinf, Isnan, Real, Rint, Round, Scale, Sign, Trunc};
63pub use pow::{Exponent, Pow, RealPow, Root};
64pub use slice::{Arg, Binary, CHUNK, Flipped, KernelItem, Results, Ternary, Unary, Wrapped, pairs};
65pub use special::{BetaincAt, Erfc, Erfcinv, Erfinv};
66pub use sqrts::{RealSqrt, RealSqrt1mx2, RealSqrtm1mx2, RealSqrtx2m1, Sqrt, Sqrtm1mx2};
67pub use trig::{Cos, Sin, Tan};
68/// The crates whose types compact buffers hold, which the exported width macros name.
69pub use {half, num_complex};
70
71/// How many items a function handles between checks for an interrupt ([`interrupt::check`]), which keeps the check out of the loop over items.
72pub const BETWEEN_CHECKS: usize = 1024;
73
74/// 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.
75macro_rules! functions {
76 ($($(#[doc = $doc:literal])* $name:ident($($x:ident),+) in $module:ident $(as $object:ident)?;)+) => {
77 $(
78 $(#[doc = $doc])*
79 #[inline(always)]
80 pub fn $name<F: Float>($($x: F),+) -> F { functions!(@call FUSED, $module, $name, ($($x),+)) }
81 )+
82 $(functions!(@object $($object)?; $name, $module, ($($x),+));)+
83 /// Each function's doc by its name, as `function_doc!(ln)`, with its lines joined by newlines.
84 #[macro_export]
85 macro_rules! function_doc {
86 $(($name) => { concat!($($doc, "\n"),*) };)+
87 }
88 /// 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.
89 pub mod fused {
90 use crate::Float;
91 $(
92 #[doc = concat!("[`", stringify!($name), "`](crate::", stringify!($name), "), with `mla` fused.")]
93 #[inline(always)]
94 pub fn $name<F: Float>($($x: F),+) -> F { functions!(@call true, $module, $name, ($($x),+)) }
95 )+
96 }
97 /// rustymath's functions with `mla` unfused, whatever the target's `FUSED`.
98 pub mod unfused {
99 use crate::Float;
100 $(
101 #[doc = concat!("[`", stringify!($name), "`](crate::", stringify!($name), "), with `mla` unfused.")]
102 #[inline(always)]
103 pub fn $name<F: Float>($($x: F),+) -> F { functions!(@call false, $module, $name, ($($x),+)) }
104 )+
105 }
106 };
107 (@object; $name:ident, $module:ident, ($($x:ident),+)) => {};
108 (@object $object:ident; $name:ident, $module:ident, ($($x:ident),+)) => {
109 #[doc = concat!("[`", stringify!($name), "`](crate::", stringify!($name), ") as a function object.")]
110 #[derive(Clone, Copy, Debug)]
111 pub struct $object;
112 };
113 (@call $fused:expr, $module:ident, $name:ident, ($($x:ident),+)) => { F::from_working(crate::$module::$name::<_, { $fused }>($($x.working()),+)) };
114}
115
116functions! {
117 /// e^x. `exp(∞)` is ∞, `exp(-∞)` is 0, and `exp(NaN)` is NaN.
118 exp(x) in exp;
119 /// e^x - 1, without the loss of digits near 0 that `exp(x) - 1` has.
120 expm1(x) in exp;
121 /// The natural logarithm of x. `ln(0)` is -∞, `ln(∞)` is ∞, and a negative real x or NaN gives NaN.
122 ln(x) in ln;
123 /// ln(1 + x), without the loss of digits near 0 that `ln(1 + x)` has. `log1p(-1)` is -∞.
124 log1p(x) in ln;
125 /// The logarithm of x to base 2. It is k exactly at 2^k.
126 log2(x) in ln;
127 /// 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.
128 log10(x) in ln;
129 /// 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.
130 logb(b, x) in ln;
131 /// 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.
132 sin(x) in trig;
133 /// 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.
134 cos(x) in trig;
135 /// 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.
136 tan(x) in trig;
137 /// asin x. A real x gives a result from -π/2 to π/2 and NaN beyond a magnitude of 1.
138 asin(x) in arc as Asin;
139 /// acos x. A real x gives a result from 0 to π and NaN beyond a magnitude of 1.
140 acos(x) in arc as Acos;
141 /// atan x. A real x gives a result from -π/2 to π/2.
142 atan(x) in arc as Atan;
143 /// 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.
144 atan2(y, x) in arc as Atan2;
145 /// sinh x. It stays finite up to where sinh itself overflows, beyond the point where e^x does.
146 sinh(x) in hyperbolic as Sinh;
147 /// cosh x. It stays finite up to where cosh itself overflows, beyond the point where e^x does.
148 cosh(x) in hyperbolic as Cosh;
149 /// tanh x.
150 tanh(x) in hyperbolic as Tanh;
151 /// asinh x.
152 asinh(x) in hyperbolic as Asinh;
153 /// acosh x. A real x below 1 gives NaN.
154 acosh(x) in hyperbolic as Acosh;
155 /// atanh x. ±1 gives ±∞, and a real x of magnitude beyond 1 gives NaN.
156 atanh(x) in hyperbolic as Atanh;
157 /// 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.
158 pow(x, y) in pow;
159 /// 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.
160 root(n, y) in pow;
161 /// The cube root of x.
162 cbrt(x) in pow as Cbrt;
163 /// √(x² + y²), with no overflow or underflow in between. An infinite argument gives ∞, even beside NaN.
164 hypot(x, y) in pow as Hypot;
165 /// √(1 - x²), without the digits that squaring x first loses near |x| = 1. A real x of magnitude beyond 1 gives NaN.
166 sqrt1mx2(x) in sqrts as Sqrt1mx2;
167 /// √(1 + x²), with no overflow where x² would overflow.
168 sqrt1px2(x) in sqrts as Sqrt1px2;
169 /// √(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.
170 sqrtx2m1(x) in sqrts as Sqrtx2m1;
171 /// Γ(x). Γ(±0) is ±∞, a negative whole x and -∞ give NaN, and Γ(∞) is ∞.
172 gamma(x) in special as Gamma;
173 /// Γ(x + 1), which is x! for a whole x from 0.
174 factorial(x) in special as Factorial;
175 /// 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.
176 lgamma(x) in special as Lgamma;
177 /// ψ(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).
178 digamma(x) in special as Digamma;
179 /// ψ⁽ⁿ⁾(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.
180 polygamma(n, x) in special as Polygamma;
181 /// 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.
182 beta(a, b) in special as Beta;
183 /// 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.
184 binom(n, k) in special as Binom;
185 /// ∂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.
186 binom_dn(n, k) in special as BinomDn;
187 /// ∂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.
188 binom_dk(n, k) in special as BinomDk;
189 /// The error function, 2/√π times the integral of e^-t² from 0 to x.
190 erf(x) in special as Erf;
191 /// 1 - erf(x), without the loss of digits for large x that `1 - erf(x)` has.
192 erfc(x) in special;
193 /// The x with erf(x) = p. ±1 gives ±∞, and a magnitude beyond 1 gives NaN.
194 erfinv(p) in special;
195 /// The x with erfc(x) = q. 0 gives ∞, 2 gives -∞, and q outside 0 to 2 gives NaN.
196 erfcinv(q) in special;
197 /// 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.
198 lambertw(x) in lambert as LambertW;
199 /// 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.
200 gammainc(a, x) in special as Gammainc;
201 /// 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.
202 gammaincc(a, x) in special as Gammaincc;
203 /// The x with P(a, x) = p, for p from 0 to 1.
204 gammaincinv(a, p) in special as Gammaincinv;
205 /// 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.
206 betainc(a, b, x) in special as Betainc;
207 /// The x with I_x(a, b) = p, for p from 0 to 1.
208 betaincinv(a, b, p) in special as Betaincinv;
209}