Source file src/pkg/math/cbrt.go
1 // Copyright 2009 The Go Authors. All rights reserved. 2 // Use of this source code is governed by a BSD-style 3 // license that can be found in the LICENSE file. 4 5 package math 6 7 /* 8 The algorithm is based in part on "Optimal Partitioning of 9 Newton's Method for Calculating Roots", by Gunter Meinardus 10 and G. D. Taylor, Mathematics of Computation © 1980 American 11 Mathematical Society. 12 (http://www.jstor.org/stable/2006387?seq=9, accessed 11-Feb-2010) 13 */ 14 15 // Cbrt returns the cube root of its argument. 16 // 17 // Special cases are: 18 // Cbrt(±0) = ±0 19 // Cbrt(±Inf) = ±Inf 20 // Cbrt(NaN) = NaN 21 func Cbrt(x float64) float64 { 22 const ( 23 A1 = 1.662848358e-01 24 A2 = 1.096040958e+00 25 A3 = 4.105032829e-01 26 A4 = 5.649335816e-01 27 B1 = 2.639607233e-01 28 B2 = 8.699282849e-01 29 B3 = 1.629083358e-01 30 B4 = 2.824667908e-01 31 C1 = 4.190115298e-01 32 C2 = 6.904625373e-01 33 C3 = 6.46502159e-02 34 C4 = 1.412333954e-01 35 ) 36 // special cases 37 switch { 38 case x == 0 || IsNaN(x) || IsInf(x, 0): 39 return x 40 } 41 sign := false 42 if x < 0 { 43 x = -x 44 sign = true 45 } 46 // Reduce argument and estimate cube root 47 f, e := Frexp(x) // 0.5 <= f < 1.0 48 m := e % 3 49 if m > 0 { 50 m -= 3 51 e -= m // e is multiple of 3 52 } 53 switch m { 54 case 0: // 0.5 <= f < 1.0 55 f = A1*f + A2 - A3/(A4+f) 56 case -1: 57 f *= 0.5 // 0.25 <= f < 0.5 58 f = B1*f + B2 - B3/(B4+f) 59 default: // m == -2 60 f *= 0.25 // 0.125 <= f < 0.25 61 f = C1*f + C2 - C3/(C4+f) 62 } 63 y := Ldexp(f, e/3) // e/3 = exponent of cube root 64 65 // Iterate 66 s := y * y * y 67 t := s + x 68 y *= (t + x) / (s + t) 69 // Reiterate 70 s = (y*y*y - x) / x 71 y -= y * (((14.0/81.0)*s-(2.0/9.0))*s + (1.0 / 3.0)) * s 72 if sign { 73 y = -y 74 } 75 return y 76 }