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         return libc_pow(x, y)
  40 }
  41
  42 //extern pow
  43 func libc_pow(float64, float64) float64
  44
  45 func pow(x, y float64) float64 {
  46         switch {
  47         case y == 0 || x == 1:
  48                 return 1
  49         case y == 1:
  50                 return x
  51         case IsNaN(x) || IsNaN(y):
  52                 return NaN()
  53         case x == 0:
  54                 switch {
  55                 case y < 0:
  56                         if isOddInt(y) {
  57                                 return Copysign(Inf(1), x)
  58                         }
  59                         return Inf(1)
  60                 case y > 0:
  61                         if isOddInt(y) {
  62                                 return x
  63                         }
  64                         return 0
  65                 }
  66         case IsInf(y, 0):
  67                 switch {
  68                 case x == -1:
  69                         return 1
  70                 case (Abs(x) < 1) == IsInf(y, 1):
  71                         return 0
  72                 default:
  73                         return Inf(1)
  74                 }
  75         case IsInf(x, 0):
  76                 if IsInf(x, -1) {
  77                         return Pow(1/x, -y) // Pow(-0, -y)
  78                 }
  79                 switch {
  80                 case y < 0:
  81                         return 0
  82                 case y > 0:
  83                         return Inf(1)
  84                 }
  85         case y == 0.5:
  86                 return Sqrt(x)
  87         case y == -0.5:
  88                 return 1 / Sqrt(x)
  89         }
  90
  91         yi, yf := Modf(Abs(y))
  92         if yf != 0 && x < 0 {
  93                 return NaN()
  94         }
  95         if yi >= 1<<63 {
  96                 // yi is a large even int that will lead to overflow (or underflow to 0)
  97                 // for all x except -1 (x == 1 was handled earlier)
  98                 switch {
  99                 case x == -1:
 100                         return 1
 101                 case (Abs(x) < 1) == (y > 0):
 102                         return 0
 103                 default:
 104                         return Inf(1)
 105                 }
 106         }
 107
 108         // ans = a1 * 2**ae (= 1 for now).
 109         a1 := 1.0
 110         ae := 0
 111
 112         // ans *= x**yf
 113         if yf != 0 {
 114                 if yf > 0.5 {
 115                         yf--
 116                         yi++
 117                 }
 118                 a1 = Exp(yf * Log(x))
 119         }
 120
 121         // ans *= x**yi
 122         // by multiplying in successive squarings
 123         // of x according to bits of yi.
 124         // accumulate powers of two into exp.
 125         x1, xe := Frexp(x)
 126         for i := int64(yi); i != 0; i >>= 1 {
 127                 if xe < -1<<12 || 1<<12 < xe {
 128                         // catch xe before it overflows the left shift below
 129                         // Since i !=0 it has at least one bit still set, so ae will accumulate xe
 130                         // on at least one more iteration, ae += xe is a lower bound on ae
 131                         // the lower bound on ae exceeds the size of a float64 exp
 132                         // so the final call to Ldexp will produce under/overflow (0/Inf)
 133                         ae += xe
 134                         break
 135                 }
 136                 if i&1 == 1 {
 137                         a1 *= x1
 138                         ae += xe
 139                 }
 140                 x1 *= x1
 141                 xe <<= 1
 142                 if x1 < .5 {
 143                         x1 += x1
 144                         xe--
 145                 }
 146         }
 147
 148         // ans = a1*2**ae
 149         // if y < 0 { ans = 1 / ans }
 150         // but in the opposite order
 151         if y < 0 {
 152                 a1 = 1 / a1
 153                 ae = -ae
 154         }
 155         return Ldexp(a1, ae)
 156 }