deep_time/math/
k_tan.rs

1// origin: FreeBSD /usr/src/lib/msun/src/k_tan.c
2//
3// ====================================================
4// Copyright 2004 Sun Microsystems, Inc.  All Rights Reserved.
5//
6// Permission to use, copy, modify, and distribute this
7// software is freely granted, provided that this notice
8// is preserved.
9// ====================================================
10
11#![allow(clippy::indexing_slicing)]
12#![allow(clippy::excessive_precision)]
13#![allow(clippy::approx_constant)]
14#![allow(clippy::eq_op)]
15
16use crate::Real;
17
18static T: [Real; 13] = [
19    Real::from_bits(0x3fd5555555555563),
20    Real::from_bits(0x3fc111111110fe7a),
21    Real::from_bits(0x3faba1ba1bb341fe),
22    Real::from_bits(0x3f9664f48406d637),
23    Real::from_bits(0x3f8226e3e96e8493),
24    Real::from_bits(0x3f6d6d22c9560328),
25    Real::from_bits(0x3f57dbc8fee08315),
26    Real::from_bits(0x3f4344d8f2f26501),
27    Real::from_bits(0x3f3026f71a8d1068),
28    Real::from_bits(0x3f147e88a03792a6),
29    Real::from_bits(0x3f12b80f32f0a7e9),
30    Real::from_bits(0xbef375cbdb605373),
31    Real::from_bits(0x3efb2a7074bf7ad4),
32];
33
34const ONE: Real = Real::from_bits(0x3ff0000000000000);
35const TWO: Real = Real::from_bits(0x4000000000000000);
36const PIO4: Real = Real::from_bits(0x3fe921fb54442d18);
37const PIO4_LO: Real = Real::from_bits(0x3c81a62633145c07);
38
39/// Exact port of FreeBSD msun `k_tan` using the **original 0/1 calling convention**
40/// that `s_tan.c` (and the top-level `tan`) expects.
41///
42/// `odd == 0` → return tan(x+y)
43/// `odd == 1` → return -1/tan(x+y)
44pub(crate) const fn k_tan(mut x: Real, mut y: Real, odd: i32) -> Real {
45    let hx = (Real::to_bits(x) >> 32) as u32;
46    let big = (hx & 0x7fffffff) >= 0x3fe59428; /* |x| >= 0.6744 */
47
48    if big {
49        let sign = hx >> 31;
50        if sign != 0 {
51            x = -x;
52            y = -y;
53        }
54        x = (PIO4 - x) + (PIO4_LO - y);
55        y = Real::from_bits(0);
56    }
57
58    let z = x * x;
59    let w = z * z;
60
61    let r = T[1] + w * (T[3] + w * (T[5] + w * (T[7] + w * (T[9] + w * T[11]))));
62    let v = z * (T[2] + w * (T[4] + w * (T[6] + w * (T[8] + w * (T[10] + w * T[12])))));
63
64    let s = z * x;
65    let r = y + z * (s * (r + v) + y) + s * T[0];
66    let w = x + r;
67
68    if big {
69        let sign = hx >> 31; // ORIGINAL sign (important!)
70
71        // This is the exact line from the original FreeBSD source.
72        // We avoid any From<i64> by using a simple if (odd is only 0 or 1).
73        let s_val = if odd == 0 { ONE } else { ONE - TWO };
74
75        let tmp = r - (w * w / (w + s_val));
76        let v = s_val - TWO * (x + tmp);
77
78        return if sign != 0 { -v } else { v };
79    }
80
81    if odd == 0 {
82        return w;
83    }
84
85    /* -1.0/(x+r) with extra precision */
86    let w0 = zero_low_word(w);
87    let v = r - (w0 - x);
88    let a = -ONE / w;
89    let a0 = zero_low_word(a);
90    a0 + a * (ONE + a0 * w0 + a0 * v)
91}
92
93#[inline]
94const fn zero_low_word(x: Real) -> Real {
95    Real::from_bits(Real::to_bits(x) & 0xffff_ffff_0000_0000)
96}
deep_time/math/k_tan.rs

deep_time/math/
k_tan.rs
