github.com/ccccaoqing/test@v0.0.0-20220510085219-3985d23445c0/src/math/cbrt.go (about)

     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 x.
    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  }