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  // The go code is a modified version of the original C code from
     8  // http://www.netlib.org/fdlibm/s_cbrt.c and came with this notice.
     9  //
    10  // ====================================================
    11  // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
    12  //
    13  // Developed at SunSoft, a Sun Microsystems, Inc. business.
    14  // Permission to use, copy, modify, and distribute this
    15  // software is freely granted, provided that this notice
    16  // is preserved.
    17  // ====================================================
    18  
    19  // Cbrt returns the cube root of x.
    20  //
    21  // Special cases are:
    22  //
    23  //	Cbrt(±0) = ±0
    24  //	Cbrt(±Inf) = ±Inf
    25  //	Cbrt(NaN) = NaN
    26  func Cbrt(x float64) float64 {
    27  	if haveArchCbrt {
    28  		return archCbrt(x)
    29  	}
    30  	return cbrt(x)
    31  }
    32  
    33  func cbrt(x float64) float64 {
    34  	const (
    35  		B1             = 715094163                   // (682-0.03306235651)*2**20
    36  		B2             = 696219795                   // (664-0.03306235651)*2**20
    37  		C              = 5.42857142857142815906e-01  // 19/35     = 0x3FE15F15F15F15F1
    38  		D              = -7.05306122448979611050e-01 // -864/1225 = 0xBFE691DE2532C834
    39  		E              = 1.41428571428571436819e+00  // 99/70     = 0x3FF6A0EA0EA0EA0F
    40  		F              = 1.60714285714285720630e+00  // 45/28     = 0x3FF9B6DB6DB6DB6E
    41  		G              = 3.57142857142857150787e-01  // 5/14      = 0x3FD6DB6DB6DB6DB7
    42  		SmallestNormal = 2.22507385850720138309e-308 // 2**-1022  = 0x0010000000000000
    43  	)
    44  	// special cases
    45  	switch {
    46  	case x == 0 || IsNaN(x) || IsInf(x, 0):
    47  		return x
    48  	}
    49  
    50  	sign := false
    51  	if x < 0 {
    52  		x = -x
    53  		sign = true
    54  	}
    55  
    56  	// rough cbrt to 5 bits
    57  	t := Float64frombits(Float64bits(x)/3 + B1<<32)
    58  	if x < SmallestNormal {
    59  		// subnormal number
    60  		t = float64(1 << 54) // set t= 2**54
    61  		t *= x
    62  		t = Float64frombits(Float64bits(t)/3 + B2<<32)
    63  	}
    64  
    65  	// new cbrt to 23 bits
    66  	r := t * t / x
    67  	s := C + r*t
    68  	t *= G + F/(s+E+D/s)
    69  
    70  	// chop to 22 bits, make larger than cbrt(x)
    71  	t = Float64frombits(Float64bits(t)&(0xFFFFFFFFC<<28) + 1<<30)
    72  
    73  	// one step newton iteration to 53 bits with error less than 0.667ulps
    74  	s = t * t // t*t is exact
    75  	r = x / s
    76  	w := t + t
    77  	r = (r - t) / (w + r) // r-s is exact
    78  	t = t + t*r
    79  
    80  	// restore the sign bit
    81  	if sign {
    82  		t = -t
    83  	}
    84  	return t
    85  }
    86
View as plain text