v

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
}
65