libgo/go/math/pow.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 func isOddInt(x float64) bool {
   8         xi, xf := Modf(x)
   9         return xf == 0 && int64(xi)&1 == 1
  10 }
  11
  12 // Special cases taken from FreeBSD's /usr/src/lib/msun/src/e_pow.c
  13 // updated by IEEE Std. 754-2008 "Section 9.2.1 Special values".
  14
  15 // Pow returns x**y, the base-x exponential of y.
  16 //
  17 // Special cases are (in order):
  18 //      Pow(x, ±0) = 1 for any x
  19 //      Pow(1, y) = 1 for any y
  20 //      Pow(x, 1) = x for any x
  21 //      Pow(NaN, y) = NaN
  22 //      Pow(x, NaN) = NaN
  23 //      Pow(±0, y) = ±Inf for y an odd integer < 0
  24 //      Pow(±0, -Inf) = +Inf
  25 //      Pow(±0, +Inf) = +0
  26 //      Pow(±0, y) = +Inf for finite y < 0 and not an odd integer
  27 //      Pow(±0, y) = ±0 for y an odd integer > 0
  28 //      Pow(±0, y) = +0 for finite y > 0 and not an odd integer
  29 //      Pow(-1, ±Inf) = 1
  30 //      Pow(x, +Inf) = +Inf for |x| > 1
  31 //      Pow(x, -Inf) = +0 for |x| > 1
  32 //      Pow(x, +Inf) = +0 for |x| < 1
  33 //      Pow(x, -Inf) = +Inf for |x| < 1
  34 //      Pow(+Inf, y) = +Inf for y > 0
  35 //      Pow(+Inf, y) = +0 for y < 0
  36 //      Pow(-Inf, y) = Pow(-0, -y)
  37 //      Pow(x, y) = NaN for finite x < 0 and finite non-integer y
  38 func Pow(x, y float64) float64 {
  39         switch {
  40         case y == 0 || x == 1:
  41                 return 1
  42         case y == 1:
  43                 return x
  44         case y == 0.5:
  45                 return Sqrt(x)
  46         case y == -0.5:
  47                 return 1 / Sqrt(x)
  48         case IsNaN(x) || IsNaN(y):
  49                 return NaN()
  50         case x == 0:
  51                 switch {
  52                 case y < 0:
  53                         if isOddInt(y) {
  54                                 return Copysign(Inf(1), x)
  55                         }
  56                         return Inf(1)
  57                 case y > 0:
  58                         if isOddInt(y) {
  59                                 return x
  60                         }
  61                         return 0
  62                 }
  63         case IsInf(y, 0):
  64                 switch {
  65                 case x == -1:
  66                         return 1
  67                 case (Abs(x) < 1) == IsInf(y, 1):
  68                         return 0
  69                 default:
  70                         return Inf(1)
  71                 }
  72         case IsInf(x, 0):
  73                 if IsInf(x, -1) {
  74                         return Pow(1/x, -y) // Pow(-0, -y)
  75                 }
  76                 switch {
  77                 case y < 0:
  78                         return 0
  79                 case y > 0:
  80                         return Inf(1)
  81                 }
  82         }
  83
  84         absy := y
  85         flip := false
  86         if absy < 0 {
  87                 absy = -absy
  88                 flip = true
  89         }
  90         yi, yf := Modf(absy)
  91         if yf != 0 && x < 0 {
  92                 return NaN()
  93         }
  94         if yi >= 1<<63 {
  95                 return Exp(y * Log(x))
  96         }
  97
  98         // ans = a1 * 2**ae (= 1 for now).
  99         a1 := 1.0
 100         ae := 0
 101
 102         // ans *= x**yf
 103         if yf != 0 {
 104                 if yf > 0.5 {
 105                         yf--
 106                         yi++
 107                 }
 108                 a1 = Exp(yf * Log(x))
 109         }
 110
 111         // ans *= x**yi
 112         // by multiplying in successive squarings
 113         // of x according to bits of yi.
 114         // accumulate powers of two into exp.
 115         x1, xe := Frexp(x)
 116         for i := int64(yi); i != 0; i >>= 1 {
 117                 if i&1 == 1 {
 118                         a1 *= x1
 119                         ae += xe
 120                 }
 121                 x1 *= x1
 122                 xe <<= 1
 123                 if x1 < .5 {
 124                         x1 += x1
 125                         xe--
 126                 }
 127         }
 128
 129         // ans = a1*2**ae
 130         // if flip { ans = 1 / ans }
 131         // but in the opposite order
 132         if flip {
 133                 a1 = 1 / a1
 134                 ae = -ae
 135         }
 136         return Ldexp(a1, ae)
 137 }