libm/math/
sincos.rs

Help
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#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
16pub fn sincos(x: f64) -> (f64, f64) {
17    let s: f64;
18    let c: f64;
19    let mut ix: u32;
20
21    ix = get_high_word(x);
22    ix &= 0x7fffffff;
23
24    /* |x| ~< pi/4 */
25    if ix <= 0x3fe921fb {
26        /* if |x| < 2**-27 * sqrt(2) */
27        if ix < 0x3e46a09e {
28            /* raise inexact if x!=0 and underflow if subnormal */
29            let x1p120 = f64::from_bits(0x4770000000000000); // 0x1p120 == 2^120
30            if ix < 0x00100000 {
31                force_eval!(x / x1p120);
32            } else {
33                force_eval!(x + x1p120);
34            }
35            return (x, 1.0);
36        }
37        return (k_sin(x, 0.0, 0), k_cos(x, 0.0));
38    }
39
40    /* sincos(Inf or NaN) is NaN */
41    if ix >= 0x7ff00000 {
42        let rv = x - x;
43        return (rv, rv);
44    }
45
46    /* argument reduction needed */
47    let (n, y0, y1) = rem_pio2(x);
48    s = k_sin(y0, y1, 1);
49    c = k_cos(y0, y1);
50    match n & 3 {
51        0 => (s, c),
52        1 => (c, -s),
53        2 => (-s, -c),
54        3 => (-c, s),
55        #[cfg(debug_assertions)]
56        _ => unreachable!(),
57        #[cfg(not(debug_assertions))]
58        _ => (0.0, 1.0),
59    }
60}
61
62// These tests are based on those from sincosf.rs
63#[cfg(test)]
64mod tests {
65    use super::sincos;
66
67    const TOLERANCE: f64 = 1e-6;
68
69    #[test]
70    fn with_pi() {
71        let (s, c) = sincos(core::f64::consts::PI);
72        assert!(
73            (s - 0.0).abs() < TOLERANCE,
74            "|{} - {}| = {} >= {}",
75            s,
76            0.0,
77            (s - 0.0).abs(),
78            TOLERANCE
79        );
80        assert!(
81            (c + 1.0).abs() < TOLERANCE,
82            "|{} + {}| = {} >= {}",
83            c,
84            1.0,
85            (s + 1.0).abs(),
86            TOLERANCE
87        );
88    }
89
90    #[test]
91    fn rotational_symmetry() {
92        use core::f64::consts::PI;
93        const N: usize = 24;
94        for n in 0..N {
95            let theta = 2. * PI * (n as f64) / (N as f64);
96            let (s, c) = sincos(theta);
97            let (s_plus, c_plus) = sincos(theta + 2. * PI);
98            let (s_minus, c_minus) = sincos(theta - 2. * PI);
99
100            assert!(
101                (s - s_plus).abs() < TOLERANCE,
102                "|{} - {}| = {} >= {}",
103                s,
104                s_plus,
105                (s - s_plus).abs(),
106                TOLERANCE
107            );
108            assert!(
109                (s - s_minus).abs() < TOLERANCE,
110                "|{} - {}| = {} >= {}",
111                s,
112                s_minus,
113                (s - s_minus).abs(),
114                TOLERANCE
115            );
116            assert!(
117                (c - c_plus).abs() < TOLERANCE,
118                "|{} - {}| = {} >= {}",
119                c,
120                c_plus,
121                (c - c_plus).abs(),
122                TOLERANCE
123            );
124            assert!(
125                (c - c_minus).abs() < TOLERANCE,
126                "|{} - {}| = {} >= {}",
127                c,
128                c_minus,
129                (c - c_minus).abs(),
130                TOLERANCE
131            );
132        }
133    }
134}
libm/math/sincos.rs

libm/math/
sincos.rs
