libgo/go/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 /*
   8         The algorithm is based in part on "Optimal Partitioning of
   9         Newton's Method for Calculating Roots", by Gunter Meinardus
  10         and G. D. Taylor, Mathematics of Computation © 1980 American
  11         Mathematical Society.
  12         (http://www.jstor.org/stable/2006387?seq=9, accessed 11-Feb-2010)
  13 */
  14
  15 // Cbrt returns the cube root of its argument.
  16 //
  17 // Special cases are:
  18 //      Cbrt(±0) = ±0
  19 //      Cbrt(±Inf) = ±Inf
  20 //      Cbrt(NaN) = NaN
  21 func Cbrt(x float64) float64 {
  22         const (
  23                 A1 = 1.662848358e-01
  24                 A2 = 1.096040958e+00
  25                 A3 = 4.105032829e-01
  26                 A4 = 5.649335816e-01
  27                 B1 = 2.639607233e-01
  28                 B2 = 8.699282849e-01
  29                 B3 = 1.629083358e-01
  30                 B4 = 2.824667908e-01
  31                 C1 = 4.190115298e-01
  32                 C2 = 6.904625373e-01
  33                 C3 = 6.46502159e-02
  34                 C4 = 1.412333954e-01
  35         )
  36         // TODO(rsc): Remove manual inlining of IsNaN, IsInf
  37         // when compiler does it for us
  38         // special cases
  39         switch {
  40         case x == 0 || x != x || x < -MaxFloat64 || x > MaxFloat64: // x == 0 || IsNaN(x) || IsInf(x, 0):
  41                 return x
  42         }
  43         sign := false
  44         if x < 0 {
  45                 x = -x
  46                 sign = true
  47         }
  48         // Reduce argument and estimate cube root
  49         f, e := Frexp(x) // 0.5 <= f < 1.0
  50         m := e % 3
  51         if m > 0 {
  52                 m -= 3
  53                 e -= m // e is multiple of 3
  54         }
  55         switch m {
  56         case 0: // 0.5 <= f < 1.0
  57                 f = A1*f + A2 - A3/(A4+f)
  58         case -1:
  59                 f *= 0.5 // 0.25 <= f < 0.5
  60                 f = B1*f + B2 - B3/(B4+f)
  61         default: // m == -2
  62                 f *= 0.25 // 0.125 <= f < 0.25
  63                 f = C1*f + C2 - C3/(C4+f)
  64         }
  65         y := Ldexp(f, e/3) // e/3 = exponent of cube root
  66
  67         // Iterate
  68         s := y * y * y
  69         t := s + x
  70         y *= (t + x) / (s + t)
  71         // Reiterate
  72         s = (y*y*y - x) / x
  73         y -= y * (((14.0/81.0)*s-(2.0/9.0))*s + (1.0 / 3.0)) * s
  74         if sign {
  75                 y = -y
  76         }
  77         return y
  78 }