rustymath_core/
arc.rs

1//! The inverse circle functions. `asin` and `acos` reduce as SLEEF's `xasin_u35` and `xacos_u35` do (see `THIRD-PARTY.md`). `atan` reduces to three intervals, as Cephes' `atan` does. The polynomials are minimax fits.
2use crate::{Atan, float::{Native, horner}};
3
4/// The coefficients of (asin s - s)/s³ as a polynomial in s², highest degree first, for |s| up to ½.
5const ASIN_64: [f64; 12] = [
6    0.03191217485978631, -0.01622411038464979, 0.019528172105463745, 0.006528010080512205, 0.012169203417632283, 0.01388523674904418, 0.017359704628917373,
7    0.02237175805979482, 0.030381959136550615, 0.0446428568283393, 0.07500000000337, 0.16666666666665408,
8];
9const ASIN_32: [f64; 5] = [0.04216631501913071, 0.02417951077222824, 0.04547037556767464, 0.07495297491550446, 0.16666752099990845];
10/// The coefficients of (atan u - u)/u³ as a polynomial in u², highest degree first, for |u| up to tan(π/8).
11const ATAN_64: [f64; 10] = [
12    0.02125986000857051, -0.043590968108632014, 0.05692539639740812, -0.06641121338420021, 0.07690068320282498, -0.09090782414081072, 0.1111110665366684,
13    -0.14285714195263532, 0.19999999999088985, -0.3333333333333015,
14];
15const ATAN_32: [f64; 4] = [0.08053722977638245, -0.13877679407596588, 0.19977709650993347, -0.3333294987678528];
16
17/// π/2 as the nearest float and the rest. Halving or doubling both parts gives π/4 or π exactly.
18#[inline]
19pub(crate) fn half_pi<R: Native>() -> (R, R) {
20    let (hi, lo) = if R::WIDE { (1.5707963267948966, 6.123233995736766e-17) } else { (1.5707963705062866, -4.3711390001862426e-08) };
21    (R::narrow(hi), R::narrow(lo))
22}
23
24/// asin s, for |s| up to ½, with z = s².
25#[inline]
26fn arcsine<R: Native, const FUSED: bool>(s: R, z: R) -> R {
27    let p = if R::WIDE { horner::<R, FUSED, 12>(z, ASIN_64) } else { horner::<R, FUSED, 5>(z, ASIN_32) };
28    (s * z).mla::<FUSED>(p, s)
29}
30
31/// (s, z) for `arcsine`: x itself when |x| is up to ½, and otherwise √((1 - |x|)/2), whose arcsine is (π/2 - asin |x|)/2.
32#[inline]
33fn reduce<R: Native>(x: R) -> (R, R, bool) {
34    let a = x.abs();
35    let big = a > R::narrow(0.5);
36    let z = if big { (R::one() - a) * R::narrow(0.5) } else { x * x };
37    (if big { z.sqrt() } else { x }, z, big)
38}
39
40#[inline]
41pub(crate) fn asin<R: Native, const FUSED: bool>(x: R) -> R {
42    let (s, z, big) = reduce(x);
43    let p = arcsine::<R, FUSED>(s, z);
44    let (hi, lo) = half_pi::<R>();
45    if big { (hi - (p + p - lo)).copysign(x) } else { p }
46}
47
48/// acos x: π/2 - asin x for |x| up to ½, and otherwise 2 asin s, or π less that for a negative x.
49#[inline]
50pub(crate) fn acos<R: Native, const FUSED: bool>(x: R) -> R {
51    let (s, z, big) = reduce(x);
52    let p = arcsine::<R, FUSED>(s, z);
53    let (hi, lo) = half_pi::<R>();
54    if !big { hi - (p - lo) } else if x > R::zero() { p + p } else { (hi + hi) - ((p + p) - (lo + lo)) }
55}
56
57/// atan x. |x| reduces to u with |u| up to tan(π/8): u = |x| up to tan(π/8), u = (|x| - 1)/(|x| + 1) up to tan(3π/8), where atan |x| = π/4 + atan u, and u = -1/|x| beyond, where atan |x| = π/2 + atan u. One division serves all three.
58#[inline]
59pub(crate) fn atan<R: Native, const FUSED: bool>(x: R) -> R {
60    let t = x.abs();
61    let c = R::narrow;
62    let (mid, far) = (t > c(0.41421356237309503), t > c(2.414213562373095));
63    let num = if far { -R::one() } else if mid { t - R::one() } else { t };
64    let den = if far { t } else if mid { t + R::one() } else { R::one() };
65    let u = num / den;
66    let z = u * u;
67    let p = if R::WIDE { horner::<R, FUSED, 10>(z, ATAN_64) } else { horner::<R, FUSED, 4>(z, ATAN_32) };
68    let (hi, lo) = half_pi::<R>();
69    let (hi, lo) = if far { (hi, lo) } else if mid { (hi * c(0.5), lo * c(0.5)) } else { (R::zero(), R::zero()) };
70    (hi + ((u * z).mla::<FUSED>(p, u) + lo)).copysign(x)
71}
72
73// The Cauchy distribution's cdf runs this slice form (`distribution::through`).
74forms!(Atan at real: |_, [x]| crate::slice::Total(atan::<_, { crate::FUSED }>(x)));
75
76/// [`atan2`](crate::atan2), as atan(y/x) in the point's quadrant. It takes the arctangent of the smaller magnitude over the larger, which never overflows.
77#[inline]
78pub(crate) fn atan2<R: Native, const FUSED: bool>(y: R, x: R) -> R {
79    let (ax, ay) = (x.abs(), y.abs());
80    let steep = ay > ax;
81    let (num, den) = if steep { (ax, ay) } else { (ay, ax) };
82    // The ratio is 0 at zeros and 1 at two infinities, where the quotient is NaN.
83    let q = if den == R::zero() { R::zero() } else if ax.is_infinite() & ay.is_infinite() { R::one() } else { num / den };
84    let a = atan::<R, FUSED>(q);
85    let (hi, lo) = half_pi::<R>();
86    let a = if steep { hi - (a - lo) } else { a };
87    let a = if x.is_sign_negative() { (hi + hi) - (a - (lo + lo)) } else { a };
88    a.copysign(y)
89}
90
91/// `atan2(0, x)`, the angle of the real number `x`, with no arctangent: 0, π where `x`'s sign bit is set, or NaN, bit for bit as [`atan2`] gives them.
92#[inline(always)]
93pub(crate) fn real_angle<R: Native>(x: R) -> R {
94    let (hi, lo) = half_pi::<R>();
95    if x.is_nan() { x.abs() } else if x.is_sign_negative() { (hi + hi) - (R::zero() - (lo + lo)) } else { R::zero() }
96}
97
98/// The real part of asin x, or of its principal value at a complex z, as a function object: asin x for x from -1 to 1, and asin of the nearer of -1 and 1 beyond. The real part of `1○⁻¹`.
99#[derive(Clone, Copy, Debug)]
100pub struct RealAsin;
101
102/// The real part of acos x, or of its principal value at a complex z, as a function object: acos x for x from -1 to 1, and acos of the nearer of -1 and 1 beyond. The real part of `2○⁻¹`.
103#[derive(Clone, Copy, Debug)]
104pub struct RealAcos;
rustymath_core/arc.rs

rustymath_core/
arc.rs
