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