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)] |
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 ] |
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 | } |
135 | |