/rust/registry/src/index.crates.io-1949cf8c6b5b557f/libm-0.2.16/src/math/cbrtf.rs
Line | Count | Source |
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 | | const B1: u32 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */ |
21 | | const B2: u32 = 642849266; /* B2 = (127-127.0/3-24/3-0.03306235651)*2**23 */ |
22 | | |
23 | | /// Cube root (f32) |
24 | | /// |
25 | | /// Computes the cube root of the argument. |
26 | | #[cfg_attr(assert_no_panic, no_panic::no_panic)] |
27 | 0 | pub fn cbrtf(x: f32) -> f32 { |
28 | 0 | let x1p24 = f32::from_bits(0x4b800000); // 0x1p24f === 2 ^ 24 |
29 | | |
30 | | let mut r: f64; |
31 | | let mut t: f64; |
32 | 0 | let mut ui: u32 = x.to_bits(); |
33 | 0 | let mut hx: u32 = ui & 0x7fffffff; |
34 | | |
35 | 0 | if hx >= 0x7f800000 { |
36 | | /* cbrt(NaN,INF) is itself */ |
37 | 0 | return x + x; |
38 | 0 | } |
39 | | |
40 | | /* rough cbrt to 5 bits */ |
41 | 0 | if hx < 0x00800000 { |
42 | | /* zero or subnormal? */ |
43 | 0 | if hx == 0 { |
44 | 0 | return x; /* cbrt(+-0) is itself */ |
45 | 0 | } |
46 | 0 | ui = (x * x1p24).to_bits(); |
47 | 0 | hx = ui & 0x7fffffff; |
48 | 0 | hx = hx / 3 + B2; |
49 | 0 | } else { |
50 | 0 | hx = hx / 3 + B1; |
51 | 0 | } |
52 | 0 | ui &= 0x80000000; |
53 | 0 | ui |= hx; |
54 | | |
55 | | /* |
56 | | * First step Newton iteration (solving t*t-x/t == 0) to 16 bits. In |
57 | | * double precision so that its terms can be arranged for efficiency |
58 | | * without causing overflow or underflow. |
59 | | */ |
60 | 0 | t = f32::from_bits(ui) as f64; |
61 | 0 | r = t * t * t; |
62 | 0 | t = t * (x as f64 + x as f64 + r) / (x as f64 + r + r); |
63 | | |
64 | | /* |
65 | | * Second step Newton iteration to 47 bits. In double precision for |
66 | | * efficiency and accuracy. |
67 | | */ |
68 | 0 | r = t * t * t; |
69 | 0 | t = t * (x as f64 + x as f64 + r) / (x as f64 + r + r); |
70 | | |
71 | | /* rounding to 24 bits is perfect in round-to-nearest mode */ |
72 | 0 | t as f32 |
73 | 0 | } |