1 | /* origin: FreeBSD /usr/src/lib/msun/src/e_sqrtf.c */ |
2 | /* |
3 | * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
4 | */ |
5 | /* |
6 | * ==================================================== |
7 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
8 | * |
9 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
10 | * Permission to use, copy, modify, and distribute this |
11 | * software is freely granted, provided that this notice |
12 | * is preserved. |
13 | * ==================================================== |
14 | */ |
15 | |
16 | #[cfg_attr (all(test, assert_no_panic), no_panic::no_panic)] |
17 | pub fn sqrtf(x: f32) -> f32 { |
18 | // On wasm32 we know that LLVM's intrinsic will compile to an optimized |
19 | // `f32.sqrt` native instruction, so we can leverage this for both code size |
20 | // and speed. |
21 | llvm_intrinsically_optimized! { |
22 | #[cfg(target_arch = "wasm32" )] { |
23 | return if x < 0.0 { |
24 | ::core::f32::NAN |
25 | } else { |
26 | unsafe { ::core::intrinsics::sqrtf32(x) } |
27 | } |
28 | } |
29 | } |
30 | #[cfg (target_feature = "sse" )] |
31 | { |
32 | // Note: This path is unlikely since LLVM will usually have already |
33 | // optimized sqrt calls into hardware instructions if sse is available, |
34 | // but if someone does end up here they'll apprected the speed increase. |
35 | #[cfg (target_arch = "x86" )] |
36 | use core::arch::x86::*; |
37 | #[cfg (target_arch = "x86_64" )] |
38 | use core::arch::x86_64::*; |
39 | unsafe { |
40 | let m = _mm_set_ss(x); |
41 | let m_sqrt = _mm_sqrt_ss(m); |
42 | _mm_cvtss_f32(m_sqrt) |
43 | } |
44 | } |
45 | #[cfg (not(target_feature = "sse" ))] |
46 | { |
47 | const TINY: f32 = 1.0e-30; |
48 | |
49 | let mut z: f32; |
50 | let sign: i32 = 0x80000000u32 as i32; |
51 | let mut ix: i32; |
52 | let mut s: i32; |
53 | let mut q: i32; |
54 | let mut m: i32; |
55 | let mut t: i32; |
56 | let mut i: i32; |
57 | let mut r: u32; |
58 | |
59 | ix = x.to_bits() as i32; |
60 | |
61 | /* take care of Inf and NaN */ |
62 | if (ix as u32 & 0x7f800000) == 0x7f800000 { |
63 | return x * x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */ |
64 | } |
65 | |
66 | /* take care of zero */ |
67 | if ix <= 0 { |
68 | if (ix & !sign) == 0 { |
69 | return x; /* sqrt(+-0) = +-0 */ |
70 | } |
71 | if ix < 0 { |
72 | return (x - x) / (x - x); /* sqrt(-ve) = sNaN */ |
73 | } |
74 | } |
75 | |
76 | /* normalize x */ |
77 | m = ix >> 23; |
78 | if m == 0 { |
79 | /* subnormal x */ |
80 | i = 0; |
81 | while ix & 0x00800000 == 0 { |
82 | ix <<= 1; |
83 | i = i + 1; |
84 | } |
85 | m -= i - 1; |
86 | } |
87 | m -= 127; /* unbias exponent */ |
88 | ix = (ix & 0x007fffff) | 0x00800000; |
89 | if m & 1 == 1 { |
90 | /* odd m, double x to make it even */ |
91 | ix += ix; |
92 | } |
93 | m >>= 1; /* m = [m/2] */ |
94 | |
95 | /* generate sqrt(x) bit by bit */ |
96 | ix += ix; |
97 | q = 0; |
98 | s = 0; |
99 | r = 0x01000000; /* r = moving bit from right to left */ |
100 | |
101 | while r != 0 { |
102 | t = s + r as i32; |
103 | if t <= ix { |
104 | s = t + r as i32; |
105 | ix -= t; |
106 | q += r as i32; |
107 | } |
108 | ix += ix; |
109 | r >>= 1; |
110 | } |
111 | |
112 | /* use floating add to find out rounding direction */ |
113 | if ix != 0 { |
114 | z = 1.0 - TINY; /* raise inexact flag */ |
115 | if z >= 1.0 { |
116 | z = 1.0 + TINY; |
117 | if z > 1.0 { |
118 | q += 2; |
119 | } else { |
120 | q += q & 1; |
121 | } |
122 | } |
123 | } |
124 | |
125 | ix = (q >> 1) + 0x3f000000; |
126 | ix += m << 23; |
127 | f32::from_bits(ix as u32) |
128 | } |
129 | } |
130 | |
131 | // PowerPC tests are failing on LLVM 13: https://github.com/rust-lang/rust/issues/88520 |
132 | #[cfg (not(target_arch = "powerpc64" ))] |
133 | #[cfg (test)] |
134 | mod tests { |
135 | use super::*; |
136 | use core::f32::*; |
137 | |
138 | #[test ] |
139 | fn sanity_check() { |
140 | assert_eq!(sqrtf(100.0), 10.0); |
141 | assert_eq!(sqrtf(4.0), 2.0); |
142 | } |
143 | |
144 | /// The spec: https://en.cppreference.com/w/cpp/numeric/math/sqrt |
145 | #[test ] |
146 | fn spec_tests() { |
147 | // Not Asserted: FE_INVALID exception is raised if argument is negative. |
148 | assert!(sqrtf(-1.0).is_nan()); |
149 | assert!(sqrtf(NAN).is_nan()); |
150 | for f in [0.0, -0.0, INFINITY].iter().copied() { |
151 | assert_eq!(sqrtf(f), f); |
152 | } |
153 | } |
154 | |
155 | #[test ] |
156 | fn conformance_tests() { |
157 | let values = [ |
158 | 3.14159265359f32, |
159 | 10000.0f32, |
160 | f32::from_bits(0x0000000f), |
161 | INFINITY, |
162 | ]; |
163 | let results = [1071833029u32, 1120403456u32, 456082799u32, 2139095040u32]; |
164 | |
165 | for i in 0..values.len() { |
166 | let bits = f32::to_bits(sqrtf(values[i])); |
167 | assert_eq!(results[i], bits); |
168 | } |
169 | } |
170 | } |
171 | |