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 {
71            if sign {
72                s = -k_sinf(x as f64 + S2PIO2);
73                c = -k_cosf(x as f64 + S2PIO2);
74            } else {
75                s = -k_sinf(x as f64 - S2PIO2);
76                c = -k_cosf(x as f64 - S2PIO2);
77            }
78        }
79
80        return (s, c);
81    }
82
83    /* |x| ~<= 9*pi/4 */
84    if ix <= 0x40e231d5 {
85        if ix <= 0x40afeddf {
86            /* |x| ~<= 7*pi/4 */
87            if sign {
88                s = k_cosf(x as f64 + S3PIO2);
89                c = -k_sinf(x as f64 + S3PIO2);
90            } else {
91                s = -k_cosf(x as f64 - S3PIO2);
92                c = k_sinf(x as f64 - S3PIO2);
93            }
94        } else {
95            if sign {
96                s = k_sinf(x as f64 + S4PIO2);
97                c = k_cosf(x as f64 + S4PIO2);
98            } else {
99                s = k_sinf(x as f64 - S4PIO2);
100                c = k_cosf(x as f64 - S4PIO2);
101            }
102        }
103
104        return (s, c);
105    }
106
107    /* sin(Inf or NaN) is NaN */
108    if ix >= 0x7f800000 {
109        let rv = x - x;
110        return (rv, rv);
111    }
112
113    /* general argument reduction needed */
114    let (n, y) = rem_pio2f(x);
115    s = k_sinf(y);
116    c = k_cosf(y);
117    match n & 3 {
118        0 => (s, c),
119        1 => (c, -s),
120        2 => (-s, -c),
121        3 => (-c, s),
122        #[cfg(debug_assertions)]
123        _ => unreachable!(),
124        #[cfg(not(debug_assertions))]
125        _ => (0.0, 1.0),
126    }
127}
128
129// PowerPC tests are failing on LLVM 13: https://github.com/rust-lang/rust/issues/88520
130#[cfg(not(target_arch = "powerpc64"))]
131#[cfg(test)]
132mod tests {
133    use super::sincosf;
134
135    #[test]
136    fn rotational_symmetry() {
137        use core::f32::consts::PI;
138        const N: usize = 24;
139        for n in 0..N {
140            let theta = 2. * PI * (n as f32) / (N as f32);
141            let (s, c) = sincosf(theta);
142            let (s_plus, c_plus) = sincosf(theta + 2. * PI);
143            let (s_minus, c_minus) = sincosf(theta - 2. * PI);
144
145            const TOLERANCE: f32 = 1e-6;
146            assert!(
147                (s - s_plus).abs() < TOLERANCE,
148                "|{} - {}| = {} >= {}",
149                s,
150                s_plus,
151                (s - s_plus).abs(),
152                TOLERANCE
153            );
154            assert!(
155                (s - s_minus).abs() < TOLERANCE,
156                "|{} - {}| = {} >= {}",
157                s,
158                s_minus,
159                (s - s_minus).abs(),
160                TOLERANCE
161            );
162            assert!(
163                (c - c_plus).abs() < TOLERANCE,
164                "|{} - {}| = {} >= {}",
165                c,
166                c_plus,
167                (c - c_plus).abs(),
168                TOLERANCE
169            );
170            assert!(
171                (c - c_minus).abs() < TOLERANCE,
172                "|{} - {}| = {} >= {}",
173                c,
174                c_minus,
175                (c - c_minus).abs(),
176                TOLERANCE
177            );
178        }
179    }
180}
libm/math/sincosf.rs

libm/math/
sincosf.rs
