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