/rust/registry/src/index.crates.io-1949cf8c6b5b557f/libm-0.2.11/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 | | 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 | 0 | pub fn cbrtf(x: f32) -> f32 { |
30 | 0 | let x1p24 = f32::from_bits(0x4b800000); // 0x1p24f === 2 ^ 24 |
31 | | |
32 | | let mut r: f64; |
33 | | let mut t: f64; |
34 | 0 | let mut ui: u32 = x.to_bits(); |
35 | 0 | let mut hx: u32 = ui & 0x7fffffff; |
36 | | |
37 | 0 | if hx >= 0x7f800000 { |
38 | | /* cbrt(NaN,INF) is itself */ |
39 | 0 | return x + x; |
40 | 0 | } |
41 | | |
42 | | /* rough cbrt to 5 bits */ |
43 | 0 | if hx < 0x00800000 { |
44 | | /* zero or subnormal? */ |
45 | 0 | if hx == 0 { |
46 | 0 | return x; /* cbrt(+-0) is itself */ |
47 | 0 | } |
48 | 0 | ui = (x * x1p24).to_bits(); |
49 | 0 | hx = ui & 0x7fffffff; |
50 | 0 | hx = hx / 3 + B2; |
51 | 0 | } else { |
52 | 0 | hx = hx / 3 + B1; |
53 | 0 | } |
54 | 0 | ui &= 0x80000000; |
55 | 0 | 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 | 0 | t = f32::from_bits(ui) as f64; |
63 | 0 | r = t * t * t; |
64 | 0 | 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 | 0 | r = t * t * t; |
71 | 0 | 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 | 0 | t as f32 |
75 | 0 | } |