...
Run Format

Source file src/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	// 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	//	Cbrt(±0) = ±0
    23	//	Cbrt(±Inf) = ±Inf
    24	//	Cbrt(NaN) = NaN
    25	func Cbrt(x float64) float64 {
    26		const (
    27			B1             = 715094163                   // (682-0.03306235651)*2**20
    28			B2             = 696219795                   // (664-0.03306235651)*2**20
    29			C              = 5.42857142857142815906e-01  // 19/35     = 0x3FE15F15F15F15F1
    30			D              = -7.05306122448979611050e-01 // -864/1225 = 0xBFE691DE2532C834
    31			E              = 1.41428571428571436819e+00  // 99/70     = 0x3FF6A0EA0EA0EA0F
    32			F              = 1.60714285714285720630e+00  // 45/28     = 0x3FF9B6DB6DB6DB6E
    33			G              = 3.57142857142857150787e-01  // 5/14      = 0x3FD6DB6DB6DB6DB7
    34			SmallestNormal = 2.22507385850720138309e-308 // 2**-1022  = 0x0010000000000000
    35		)
    36		// special cases
    37		switch {
    38		case x == 0 || IsNaN(x) || IsInf(x, 0):
    39			return x
    40		}
    41	
    42		sign := false
    43		if x < 0 {
    44			x = -x
    45			sign = true
    46		}
    47	
    48		// rough cbrt to 5 bits
    49		t := Float64frombits(Float64bits(x)/3 + B1<<32)
    50		if x < SmallestNormal {
    51			// subnormal number
    52			t = float64(1 << 54) // set t= 2**54
    53			t *= x
    54			t = Float64frombits(Float64bits(t)/3 + B2<<32)
    55		}
    56	
    57		// new cbrt to 23 bits
    58		r := t * t / x
    59		s := C + r*t
    60		t *= G + F/(s+E+D/s)
    61	
    62		// chop to 22 bits, make larger than cbrt(x)
    63		t = Float64frombits(Float64bits(t)&(0xFFFFFFFFC<<28) + 1<<30)
    64	
    65		// one step newton iteration to 53 bits with error less than 0.667ulps
    66		s = t * t // t*t is exact
    67		r = x / s
    68		w := t + t
    69		r = (r - t) / (w + r) // r-s is exact
    70		t = t + t*r
    71	
    72		// restore the sign bit
    73		if sign {
    74			t = -t
    75		}
    76		return t
    77	}
    78	

View as plain text