libm/math/
sincosf.rs

Help
1/* origin: FreeBSD /usr/src/lib/msun/src/s_sinf.c */
2/*
3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
4 * Optimized by Bruce D. Evans.
5 */
6/*
7 * ====================================================
8 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
9 *
10 * Developed at SunPro, a Sun Microsystems, Inc. business.
11 * Permission to use, copy, modify, and distribute this
12 * software is freely granted, provided that this notice
13 * is preserved.
14 * ====================================================
15 */
16
17use super::{k_cosf, k_sinf, rem_pio2f};
18
19/* Small multiples of pi/2 rounded to double precision. */
20const PI_2: f64 = 0.5 * 3.1415926535897931160E+00;
21const S1PIO2: f64 = 1.0 * PI_2; /* 0x3FF921FB, 0x54442D18 */
22const S2PIO2: f64 = 2.0 * PI_2; /* 0x400921FB, 0x54442D18 */
23const S3PIO2: f64 = 3.0 * PI_2; /* 0x4012D97C, 0x7F3321D2 */
24const S4PIO2: f64 = 4.0 * PI_2; /* 0x401921FB, 0x54442D18 */
25
26/// Both the sine and cosine of `x` (f32).
27///
28/// `x` is specified in radians and the return value is (sin(x), cos(x)).
29#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
30pub fn sincosf(x: f32) -> (f32, f32) {
31    let s: f32;
32    let c: f32;
33    let mut ix: u32;
34    let sign: bool;
35
36    ix = x.to_bits();
37    sign = (ix >> 31) != 0;
38    ix &= 0x7fffffff;
39
40    /* |x| ~<= pi/4 */
41    if ix <= 0x3f490fda {
42        /* |x| < 2**-12 */
43        if ix < 0x39800000 {
44            /* raise inexact if x!=0 and underflow if subnormal */
45
46            let x1p120 = f32::from_bits(0x7b800000); // 0x1p120 == 2^120
47            if ix < 0x00100000 {
48                force_eval!(x / x1p120);
49            } else {
50                force_eval!(x + x1p120);
51            }
52            return (x, 1.0);
53        }
54        return (k_sinf(x as f64), k_cosf(x as f64));
55    }
56
57    /* |x| ~<= 5*pi/4 */
58    if ix <= 0x407b53d1 {
59        if ix <= 0x4016cbe3 {
60            /* |x| ~<= 3pi/4 */
61            if sign {
62                s = -k_cosf(x as f64 + S1PIO2);
63                c = k_sinf(x as f64 + S1PIO2);
64            } else {
65                s = k_cosf(S1PIO2 - x as f64);
66                c = k_sinf(S1PIO2 - x as f64);
67            }
68        }
69        /* -sin(x+c) is not correct if x+c could be 0: -0 vs +0 */
70        else if sign {
71            s = -k_sinf(x as f64 + S2PIO2);
72            c = -k_cosf(x as f64 + S2PIO2);
73        } else {
74            s = -k_sinf(x as f64 - S2PIO2);
75            c = -k_cosf(x as f64 - S2PIO2);
76        }
77
78        return (s, c);
79    }
80
81    /* |x| ~<= 9*pi/4 */
82    if ix <= 0x40e231d5 {
83        if ix <= 0x40afeddf {
84            /* |x| ~<= 7*pi/4 */
85            if sign {
86                s = k_cosf(x as f64 + S3PIO2);
87                c = -k_sinf(x as f64 + S3PIO2);
88            } else {
89                s = -k_cosf(x as f64 - S3PIO2);
90                c = k_sinf(x as f64 - S3PIO2);
91            }
92        } else if sign {
93            s = k_sinf(x as f64 + S4PIO2);
94            c = k_cosf(x as f64 + S4PIO2);
95        } else {
96            s = k_sinf(x as f64 - S4PIO2);
97            c = k_cosf(x as f64 - S4PIO2);
98        }
99
100        return (s, c);
101    }
102
103    /* sin(Inf or NaN) is NaN */
104    if ix >= 0x7f800000 {
105        let rv = x - x;
106        return (rv, rv);
107    }
108
109    /* general argument reduction needed */
110    let (n, y) = rem_pio2f(x);
111    s = k_sinf(y);
112    c = k_cosf(y);
113    match n & 3 {
114        0 => (s, c),
115        1 => (c, -s),
116        2 => (-s, -c),
117        3 => (-c, s),
118        #[cfg(debug_assertions)]
119        _ => unreachable!(),
120        #[cfg(not(debug_assertions))]
121        _ => (0.0, 1.0),
122    }
123}
124
125// PowerPC tests are failing on LLVM 13: https://github.com/rust-lang/rust/issues/88520
126#[cfg(not(target_arch = "powerpc64"))]
127#[cfg(test)]
128mod tests {
129    use super::sincosf;
130
131    #[test]
132    fn rotational_symmetry() {
133        use core::f32::consts::PI;
134        const N: usize = 24;
135        for n in 0..N {
136            let theta = 2. * PI * (n as f32) / (N as f32);
137            let (s, c) = sincosf(theta);
138            let (s_plus, c_plus) = sincosf(theta + 2. * PI);
139            let (s_minus, c_minus) = sincosf(theta - 2. * PI);
140
141            const TOLERANCE: f32 = 1e-6;
142            assert!(
143                (s - s_plus).abs() < TOLERANCE,
144                "|{} - {}| = {} >= {}",
145                s,
146                s_plus,
147                (s - s_plus).abs(),
148                TOLERANCE
149            );
150            assert!(
151                (s - s_minus).abs() < TOLERANCE,
152                "|{} - {}| = {} >= {}",
153                s,
154                s_minus,
155                (s - s_minus).abs(),
156                TOLERANCE
157            );
158            assert!(
159                (c - c_plus).abs() < TOLERANCE,
160                "|{} - {}| = {} >= {}",
161                c,
162                c_plus,
163                (c - c_plus).abs(),
164                TOLERANCE
165            );
166            assert!(
167                (c - c_minus).abs() < TOLERANCE,
168                "|{} - {}| = {} >= {}",
169                c,
170                c_minus,
171                (c - c_minus).abs(),
172                TOLERANCE
173            );
174        }
175    }
176}
libm/math/sincosf.rs

libm/math/
sincosf.rs
