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 | |
17 | use super::{k_cosf, k_sinf, rem_pio2f}; |
18 | |
19 | /* Small multiples of pi/2 rounded to double precision. */ |
20 | const PI_2: f32 = 0.5 * 3.1415926535897931160E+00; |
21 | const S1PIO2: f32 = 1.0 * PI_2; /* 0x3FF921FB, 0x54442D18 */ |
22 | const S2PIO2: f32 = 2.0 * PI_2; /* 0x400921FB, 0x54442D18 */ |
23 | const S3PIO2: f32 = 3.0 * PI_2; /* 0x4012D97C, 0x7F3321D2 */ |
24 | const S4PIO2: f32 = 4.0 * PI_2; /* 0x401921FB, 0x54442D18 */ |
25 | |
26 | #[cfg_attr (all(test, assert_no_panic), no_panic::no_panic)] |
27 | pub fn sincosf(x: f32) -> (f32, f32) { |
28 | let s: f32; |
29 | let c: f32; |
30 | let mut ix: u32; |
31 | let sign: bool; |
32 | |
33 | ix = x.to_bits(); |
34 | sign = (ix >> 31) != 0; |
35 | ix &= 0x7fffffff; |
36 | |
37 | /* |x| ~<= pi/4 */ |
38 | if ix <= 0x3f490fda { |
39 | /* |x| < 2**-12 */ |
40 | if ix < 0x39800000 { |
41 | /* raise inexact if x!=0 and underflow if subnormal */ |
42 | |
43 | let x1p120 = f32::from_bits(0x7b800000); // 0x1p120 == 2^120 |
44 | if ix < 0x00100000 { |
45 | force_eval!(x / x1p120); |
46 | } else { |
47 | force_eval!(x + x1p120); |
48 | } |
49 | return (x, 1.0); |
50 | } |
51 | return (k_sinf(x as f64), k_cosf(x as f64)); |
52 | } |
53 | |
54 | /* |x| ~<= 5*pi/4 */ |
55 | if ix <= 0x407b53d1 { |
56 | if ix <= 0x4016cbe3 { |
57 | /* |x| ~<= 3pi/4 */ |
58 | if sign { |
59 | s = -k_cosf((x + S1PIO2) as f64); |
60 | c = k_sinf((x + S1PIO2) as f64); |
61 | } else { |
62 | s = k_cosf((S1PIO2 - x) as f64); |
63 | c = k_sinf((S1PIO2 - x) as f64); |
64 | } |
65 | } |
66 | /* -sin(x+c) is not correct if x+c could be 0: -0 vs +0 */ |
67 | else { |
68 | if sign { |
69 | s = -k_sinf((x + S2PIO2) as f64); |
70 | c = -k_cosf((x + S2PIO2) as f64); |
71 | } else { |
72 | s = -k_sinf((x - S2PIO2) as f64); |
73 | c = -k_cosf((x - S2PIO2) as f64); |
74 | } |
75 | } |
76 | |
77 | return (s, c); |
78 | } |
79 | |
80 | /* |x| ~<= 9*pi/4 */ |
81 | if ix <= 0x40e231d5 { |
82 | if ix <= 0x40afeddf { |
83 | /* |x| ~<= 7*pi/4 */ |
84 | if sign { |
85 | s = k_cosf((x + S3PIO2) as f64); |
86 | c = -k_sinf((x + S3PIO2) as f64); |
87 | } else { |
88 | s = -k_cosf((x - S3PIO2) as f64); |
89 | c = k_sinf((x - S3PIO2) as f64); |
90 | } |
91 | } else { |
92 | if sign { |
93 | s = k_sinf((x + S4PIO2) as f64); |
94 | c = k_cosf((x + S4PIO2) as f64); |
95 | } else { |
96 | s = k_sinf((x - S4PIO2) as f64); |
97 | c = k_cosf((x - S4PIO2) as f64); |
98 | } |
99 | } |
100 | |
101 | return (s, c); |
102 | } |
103 | |
104 | /* sin(Inf or NaN) is NaN */ |
105 | if ix >= 0x7f800000 { |
106 | let rv = x - x; |
107 | return (rv, rv); |
108 | } |
109 | |
110 | /* general argument reduction needed */ |
111 | let (n, y) = rem_pio2f(x); |
112 | s = k_sinf(y); |
113 | c = k_cosf(y); |
114 | match n & 3 { |
115 | 0 => (s, c), |
116 | 1 => (c, -s), |
117 | 2 => (-s, -c), |
118 | 3 => (-c, s), |
119 | #[cfg (debug_assertions)] |
120 | _ => unreachable!(), |
121 | #[cfg (not(debug_assertions))] |
122 | _ => (0.0, 1.0), |
123 | } |
124 | } |
125 | |
126 | // PowerPC tests are failing on LLVM 13: https://github.com/rust-lang/rust/issues/88520 |
127 | #[cfg (not(target_arch = "powerpc64" ))] |
128 | #[cfg (test)] |
129 | mod tests { |
130 | use super::sincosf; |
131 | use crate::_eqf; |
132 | |
133 | #[test ] |
134 | fn with_pi() { |
135 | let (s, c) = sincosf(core::f32::consts::PI); |
136 | _eqf(s.abs(), 0.0).unwrap(); |
137 | _eqf(c, -1.0).unwrap(); |
138 | } |
139 | |
140 | #[test ] |
141 | fn rotational_symmetry() { |
142 | use core::f32::consts::PI; |
143 | const N: usize = 24; |
144 | for n in 0..N { |
145 | let theta = 2. * PI * (n as f32) / (N as f32); |
146 | let (s, c) = sincosf(theta); |
147 | let (s_plus, c_plus) = sincosf(theta + 2. * PI); |
148 | let (s_minus, c_minus) = sincosf(theta - 2. * PI); |
149 | |
150 | const TOLERANCE: f32 = 1e-6; |
151 | assert!( |
152 | (s - s_plus).abs() < TOLERANCE, |
153 | "| {} - {}| = {} >= {}" , |
154 | s, |
155 | s_plus, |
156 | (s - s_plus).abs(), |
157 | TOLERANCE |
158 | ); |
159 | assert!( |
160 | (s - s_minus).abs() < TOLERANCE, |
161 | "| {} - {}| = {} >= {}" , |
162 | s, |
163 | s_minus, |
164 | (s - s_minus).abs(), |
165 | TOLERANCE |
166 | ); |
167 | assert!( |
168 | (c - c_plus).abs() < TOLERANCE, |
169 | "| {} - {}| = {} >= {}" , |
170 | c, |
171 | c_plus, |
172 | (c - c_plus).abs(), |
173 | TOLERANCE |
174 | ); |
175 | assert!( |
176 | (c - c_minus).abs() < TOLERANCE, |
177 | "| {} - {}| = {} >= {}" , |
178 | c, |
179 | c_minus, |
180 | (c - c_minus).abs(), |
181 | TOLERANCE |
182 | ); |
183 | } |
184 | } |
185 | } |
186 | |