1 | /* origin: FreeBSD /usr/src/lib/msun/src/s_cbrtf.c */ |
2 | /* |
3 | * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
4 | * Debugged and 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 | /* cbrtf(x) |
17 | * Return cube root of x |
18 | */ |
19 | |
20 | use core::f32; |
21 | |
22 | const B1: u32 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */ |
23 | const B2: u32 = 642849266; /* B2 = (127-127.0/3-24/3-0.03306235651)*2**23 */ |
24 | |
25 | /// Cube root (f32) |
26 | /// |
27 | /// Computes the cube root of the argument. |
28 | #[cfg_attr (all(test, assert_no_panic), no_panic::no_panic)] |
29 | pub fn cbrtf(x: f32) -> f32 { |
30 | let x1p24 = f32::from_bits(0x4b800000); // 0x1p24f === 2 ^ 24 |
31 | |
32 | let mut r: f64; |
33 | let mut t: f64; |
34 | let mut ui: u32 = x.to_bits(); |
35 | let mut hx: u32 = ui & 0x7fffffff; |
36 | |
37 | if hx >= 0x7f800000 { |
38 | /* cbrt(NaN,INF) is itself */ |
39 | return x + x; |
40 | } |
41 | |
42 | /* rough cbrt to 5 bits */ |
43 | if hx < 0x00800000 { |
44 | /* zero or subnormal? */ |
45 | if hx == 0 { |
46 | return x; /* cbrt(+-0) is itself */ |
47 | } |
48 | ui = (x * x1p24).to_bits(); |
49 | hx = ui & 0x7fffffff; |
50 | hx = hx / 3 + B2; |
51 | } else { |
52 | hx = hx / 3 + B1; |
53 | } |
54 | ui &= 0x80000000; |
55 | ui |= hx; |
56 | |
57 | /* |
58 | * First step Newton iteration (solving t*t-x/t == 0) to 16 bits. In |
59 | * double precision so that its terms can be arranged for efficiency |
60 | * without causing overflow or underflow. |
61 | */ |
62 | t = f32::from_bits(ui) as f64; |
63 | r = t * t * t; |
64 | t = t * (x as f64 + x as f64 + r) / (x as f64 + r + r); |
65 | |
66 | /* |
67 | * Second step Newton iteration to 47 bits. In double precision for |
68 | * efficiency and accuracy. |
69 | */ |
70 | r = t * t * t; |
71 | t = t * (x as f64 + x as f64 + r) / (x as f64 + r + r); |
72 | |
73 | /* rounding to 24 bits is perfect in round-to-nearest mode */ |
74 | t as f32 |
75 | } |
76 | |