libm/math/
sincos.rs

1/* origin: FreeBSD /usr/src/lib/msun/src/s_sin.c */
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13use super::{get_high_word, k_cos, k_sin, rem_pio2};
14
15/// Both the sine and cosine of `x` (f64).
16///
17/// `x` is specified in radians and the return value is (sin(x), cos(x)).
18#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
19pub fn sincos(x: f64) -> (f64, f64) {
20    let s: f64;
21    let c: f64;
22    let mut ix: u32;
23
24    ix = get_high_word(x);
25    ix &= 0x7fffffff;
26
27    /* |x| ~< pi/4 */
28    if ix <= 0x3fe921fb {
29        /* if |x| < 2**-27 * sqrt(2) */
30        if ix < 0x3e46a09e {
31            /* raise inexact if x!=0 and underflow if subnormal */
32            let x1p120 = f64::from_bits(0x4770000000000000); // 0x1p120 == 2^120
33            if ix < 0x00100000 {
34                force_eval!(x / x1p120);
35            } else {
36                force_eval!(x + x1p120);
37            }
38            return (x, 1.0);
39        }
40        return (k_sin(x, 0.0, 0), k_cos(x, 0.0));
41    }
42
43    /* sincos(Inf or NaN) is NaN */
44    if ix >= 0x7ff00000 {
45        let rv = x - x;
46        return (rv, rv);
47    }
48
49    /* argument reduction needed */
50    let (n, y0, y1) = rem_pio2(x);
51    s = k_sin(y0, y1, 1);
52    c = k_cos(y0, y1);
53    match n & 3 {
54        0 => (s, c),
55        1 => (c, -s),
56        2 => (-s, -c),
57        3 => (-c, s),
58        #[cfg(debug_assertions)]
59        _ => unreachable!(),
60        #[cfg(not(debug_assertions))]
61        _ => (0.0, 1.0),
62    }
63}
64
65// These tests are based on those from sincosf.rs
66#[cfg(test)]
67mod tests {
68    use super::sincos;
69
70    const TOLERANCE: f64 = 1e-6;
71
72    #[test]
73    fn with_pi() {
74        let (s, c) = sincos(core::f64::consts::PI);
75        assert!(
76            (s - 0.0).abs() < TOLERANCE,
77            "|{} - {}| = {} >= {}",
78            s,
79            0.0,
80            (s - 0.0).abs(),
81            TOLERANCE
82        );
83        assert!(
84            (c + 1.0).abs() < TOLERANCE,
85            "|{} + {}| = {} >= {}",
86            c,
87            1.0,
88            (s + 1.0).abs(),
89            TOLERANCE
90        );
91    }
92
93    #[test]
94    fn rotational_symmetry() {
95        use core::f64::consts::PI;
96        const N: usize = 24;
97        for n in 0..N {
98            let theta = 2. * PI * (n as f64) / (N as f64);
99            let (s, c) = sincos(theta);
100            let (s_plus, c_plus) = sincos(theta + 2. * PI);
101            let (s_minus, c_minus) = sincos(theta - 2. * PI);
102
103            assert!(
104                (s - s_plus).abs() < TOLERANCE,
105                "|{} - {}| = {} >= {}",
106                s,
107                s_plus,
108                (s - s_plus).abs(),
109                TOLERANCE
110            );
111            assert!(
112                (s - s_minus).abs() < TOLERANCE,
113                "|{} - {}| = {} >= {}",
114                s,
115                s_minus,
116                (s - s_minus).abs(),
117                TOLERANCE
118            );
119            assert!(
120                (c - c_plus).abs() < TOLERANCE,
121                "|{} - {}| = {} >= {}",
122                c,
123                c_plus,
124                (c - c_plus).abs(),
125                TOLERANCE
126            );
127            assert!(
128                (c - c_minus).abs() < TOLERANCE,
129                "|{} - {}| = {} >= {}",
130                c,
131                c_minus,
132                (c - c_minus).abs(),
133                TOLERANCE
134            );
135        }
136    }
137}
libm/math/sincos.rs

libm/math/
sincos.rs
