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