Home | History | Annotate | Download | only in math 
      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 }
    138