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		// TODO(rsc): Remove manual inlining of IsNaN, IsInf
    37		// when compiler does it for us
    38		// special cases
    39		switch {
    40		case x == 0 || x != x || x < -MaxFloat64 || x > MaxFloat64: // x == 0 || IsNaN(x) || IsInf(x, 0):
    41			return x
    42		}
    43		sign := false
    44		if x < 0 {
    45			x = -x
    46			sign = true
    47		}
    48		// Reduce argument
    49		f, e := Frexp(x)
    50		m := e % 3
    51		if m > 0 {
    52			m -= 3
    53			e -= m // e is multiple of 3
    54		}
    55		f = Ldexp(f, m) // 0.125 <= f < 1.0
    56	
    57		// Estimate cube root
    58		switch m {
    59		case 0: // 0.5 <= f < 1.0
    60			f = A1*f + A2 - A3/(A4+f)
    61		case -1: // 0.25 <= f < 0.5
    62			f = B1*f + B2 - B3/(B4+f)
    63		default: // 0.125 <= f < 0.25
    64			f = C1*f + C2 - C3/(C4+f)
    65		}
    66		y := Ldexp(f, e/3) // e/3 = exponent of cube root
    67	
    68		// Iterate
    69		s := y * y * y
    70		t := s + x
    71		y *= (t + x) / (s + t)
    72		// Reiterate
    73		s = (y*y*y - x) / x
    74		y -= y * (((14.0/81.0)*s-(2.0/9.0))*s + (1.0 / 3.0)) * s
    75		if sign {
    76			y = -y
    77		}
    78		return y
    79	}
The Go Programming Language

Source file src/pkg/math/cbrt.go
