1 | module math |
2 | |
3 | // The vlang code is a modified version of the original C code from |
4 | // http://www.netlib.org/fdlibm/s_cbrt.c and came with this notice. |
5 | // |
6 | // ==================================================== |
7 | // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
8 | // |
9 | // Developed at SunSoft, 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 | // cbrt returns the cube root of a. |
16 | // |
17 | // special cases are: |
18 | // cbrt(±0) = ±0 |
19 | // cbrt(±inf) = ±inf |
20 | // cbrt(nan) = nan |
21 | pub fn cbrt(a f64) f64 { |
22 | mut x := a |
23 | b1 := 715094163 // (682-0.03306235651)*2**20 |
24 | b2 := 696219795 // (664-0.03306235651)*2**20 |
25 | c := 5.42857142857142815906e-01 // 19/35 = 0x3FE15F15F15F15F1 |
26 | d := -7.05306122448979611050e-01 // -864/1225 = 0xBFE691DE2532C834 |
27 | e_ := 1.41428571428571436819e+00 // 99/70 = 0x3FF6A0EA0EA0EA0F |
28 | f := 1.60714285714285720630e+00 // 45/28 = 0x3FF9B6DB6DB6DB6E |
29 | g := 3.57142857142857150787e-01 // 5/14 = 0x3FD6DB6DB6DB6DB7 |
30 | smallest_normal := 2.22507385850720138309e-308 // 2**-1022 = 0x0010000000000000 |
31 | if x == 0.0 || is_nan(x) || is_inf(x, 0) { |
32 | return x |
33 | } |
34 | mut sign := false |
35 | if x < 0 { |
36 | x = -x |
37 | sign = true |
38 | } |
39 | // rough cbrt to 5 bits |
40 | mut t := f64_from_bits(f64_bits(x) / u64(3) + (u64(b1) << 32)) |
41 | if x < smallest_normal { |
42 | // subnormal number |
43 | t = f64(u64(1) << 54) // set t= 2**54 |
44 | t *= x |
45 | t = f64_from_bits(f64_bits(t) / u64(3) + (u64(b2) << 32)) |
46 | } |
47 | // new cbrt to 23 bits |
48 | mut r := t * t / x |
49 | mut s := c + r * t |
50 | t *= g + f / (s + e_ + d / s) |
51 | // chop to 22 bits, make larger than cbrt(x) |
52 | t = f64_from_bits(f64_bits(t) & (u64(0xffffffffc) << 28) + (u64(1) << 30)) |
53 | // one step newton iteration to 53 bits with error less than 0.667ulps |
54 | s = t * t // t*t is exact |
55 | r = x / s |
56 | w := t + t |
57 | r = (r - t) / (w + r) // r-s is exact |
58 | t = t + t * r |
59 | // restore the sign bit |
60 | if sign { |
61 | t = -t |
62 | } |
63 | return t |
64 | } |