codekingpro/portable-devtools

1// Copyright 2009 The Go Authors. All rights reserved.2// Use of this source code is governed by a BSD-style3// license that can be found in the LICENSE file.4 5package math6 7// The go code is a modified version of the original C code from8// 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 this15// software is freely granted, provided that this notice16// is preserved.17// ====================================================18 19// Cbrt returns the cube root of x.20//21// Special cases are:22//23//	Cbrt(±0) = ±024//	Cbrt(±Inf) = ±Inf25//	Cbrt(NaN) = NaN26func Cbrt(x float64) float64 {27	if haveArchCbrt {28		return archCbrt(x)29	}30	return cbrt(x)31}32 33func cbrt(x float64) float64 {34	const (35		B1             = 715094163                   // (682-0.03306235651)*2**2036		B2             = 696219795                   // (664-0.03306235651)*2**2037		C              = 5.42857142857142815906e-01  // 19/35     = 0x3FE15F15F15F15F138		D              = -7.05306122448979611050e-01 // -864/1225 = 0xBFE691DE2532C83439		E              = 1.41428571428571436819e+00  // 99/70     = 0x3FF6A0EA0EA0EA0F40		F              = 1.60714285714285720630e+00  // 45/28     = 0x3FF9B6DB6DB6DB6E41		G              = 3.57142857142857150787e-01  // 5/14      = 0x3FD6DB6DB6DB6DB742		SmallestNormal = 2.22507385850720138309e-308 // 2**-1022  = 0x001000000000000043	)44	// special cases45	switch {46	case x == 0 || IsNaN(x) || IsInf(x, 0):47		return x48	}49 50	sign := false51	if x < 0 {52		x = -x53		sign = true54	}55 56	// rough cbrt to 5 bits57	t := Float64frombits(Float64bits(x)/3 + B1<<32)58	if x < SmallestNormal {59		// subnormal number60		t = float64(1 << 54) // set t= 2**5461		t *= x62		t = Float64frombits(Float64bits(t)/3 + B2<<32)63	}64 65	// new cbrt to 23 bits66	r := t * t / x67	s := C + r*t68	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.667ulps74	s = t * t // t*t is exact75	r = x / s76	w := t + t77	r = (r - t) / (w + r) // r-s is exact78	t = t + t*r79 80	// restore the sign bit81	if sign {82		t = -t83	}84	return t85}86