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 5package 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 25func Cbrt(x float64) float64 { 26 return libc_cbrt(x) 27} 28 29//extern cbrt 30func libc_cbrt(float64) float64 31 32func 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