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 5package math 6 7func 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 38func Pow(x, y float64) float64 { 39 return libc_pow(x, y) 40} 41 42//extern pow 43func libc_pow(float64, float64) float64 44 45func 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} 157