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}
138