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
13 use super::{get_high_word, k_cos, k_sin, rem_pio2};
14
15 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
sincos(x: f64) -> (f64, f64)16 pub 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)]
64 mod tests {
65 use super::sincos;
66
67 const TOLERANCE: f64 = 1e-6;
68
69 #[test]
with_pi()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]
rotational_symmetry()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 }
135