1 /* origin: FreeBSD /usr/src/lib/msun/src/e_sqrtf.c */
2 /*
3  * Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
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)]
sqrtf(x: f32) -> f3217 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]
sanity_check()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]
spec_tests()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]
conformance_tests()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