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 x.
  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         // special cases
  37         switch {
  38         case x == 0 || IsNaN(x) || IsInf(x, 0):
  39                 return x
  40         }
  41         sign := false
  42         if x < 0 {
  43                 x = -x
  44                 sign = true
  45         }
  46         // Reduce argument and estimate cube root
  47         f, e := Frexp(x) // 0.5 <= f < 1.0
  48         m := e % 3
  49         if m > 0 {
  50                 m -= 3
  51                 e -= m // e is multiple of 3
  52         }
  53         switch m {
  54         case 0: // 0.5 <= f < 1.0
  55                 f = A1*f + A2 - A3/(A4+f)
  56         case -1:
  57                 f *= 0.5 // 0.25 <= f < 0.5
  58                 f = B1*f + B2 - B3/(B4+f)
  59         default: // m == -2
  60                 f *= 0.25 // 0.125 <= f < 0.25
  61                 f = C1*f + C2 - C3/(C4+f)
  62         }
  63         y := Ldexp(f, e/3) // e/3 = exponent of cube root
  64
  65         // Iterate
  66         s := y * y * y
  67         t := s + x
  68         y *= (t + x) / (s + t)
  69         // Reiterate
  70         s = (y*y*y - x) / x
  71         y -= y * (((14.0/81.0)*s-(2.0/9.0))*s + (1.0 / 3.0)) * s
  72         if sign {
  73                 y = -y
  74         }
  75         return y
  76 }