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