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 if Abs(x) >= (1 << 53) { 9 // 1 << 53 is the largest exact integer in the float64 format. 10 // Any number outside this range will be truncated before the decimal point and therefore will always be 11 // an even integer. 12 // Without this check and if x overflows int64 the int64(xi) conversion below may produce incorrect results 13 // on some architectures (and does so on arm64). See issue #57465. 14 return false 15 } 16 17 xi, xf := Modf(x) 18 return xf == 0 && int64(xi)&1 == 1 19} 20 21// Special cases taken from FreeBSD's /usr/src/lib/msun/src/e_pow.c 22// updated by IEEE Std. 754-2008 "Section 9.2.1 Special values". 23 24// Pow returns x**y, the base-x exponential of y. 25// 26// Special cases are (in order): 27// 28// Pow(x, ±0) = 1 for any x 29// Pow(1, y) = 1 for any y 30// Pow(x, 1) = x for any x 31// Pow(NaN, y) = NaN 32// Pow(x, NaN) = NaN 33// Pow(±0, y) = ±Inf for y an odd integer < 0 34// Pow(±0, -Inf) = +Inf 35// Pow(±0, +Inf) = +0 36// Pow(±0, y) = +Inf for finite y < 0 and not an odd integer 37// Pow(±0, y) = ±0 for y an odd integer > 0 38// Pow(±0, y) = +0 for finite y > 0 and not an odd integer 39// Pow(-1, ±Inf) = 1 40// Pow(x, +Inf) = +Inf for |x| > 1 41// Pow(x, -Inf) = +0 for |x| > 1 42// Pow(x, +Inf) = +0 for |x| < 1 43// Pow(x, -Inf) = +Inf for |x| < 1 44// Pow(+Inf, y) = +Inf for y > 0 45// Pow(+Inf, y) = +0 for y < 0 46// Pow(-Inf, y) = Pow(-0, -y) 47// Pow(x, y) = NaN for finite x < 0 and finite non-integer y 48func Pow(x, y float64) float64 { 49 if haveArchPow { 50 return archPow(x, y) 51 } 52 return pow(x, y) 53} 54 55func pow(x, y float64) float64 { 56 switch { 57 case y == 0 || x == 1: 58 return 1 59 case y == 1: 60 return x 61 case IsNaN(x) || IsNaN(y): 62 return NaN() 63 case x == 0: 64 switch { 65 case y < 0: 66 if Signbit(x) && isOddInt(y) { 67 return Inf(-1) 68 } 69 return Inf(1) 70 case y > 0: 71 if Signbit(x) && isOddInt(y) { 72 return x 73 } 74 return 0 75 } 76 case IsInf(y, 0): 77 switch { 78 case x == -1: 79 return 1 80 case (Abs(x) < 1) == IsInf(y, 1): 81 return 0 82 default: 83 return Inf(1) 84 } 85 case IsInf(x, 0): 86 if IsInf(x, -1) { 87 return Pow(1/x, -y) // Pow(-0, -y) 88 } 89 switch { 90 case y < 0: 91 return 0 92 case y > 0: 93 return Inf(1) 94 } 95 case y == 0.5: 96 return Sqrt(x) 97 case y == -0.5: 98 return 1 / Sqrt(x) 99 } 100 101 yi, yf := Modf(Abs(y)) 102 if yf != 0 && x < 0 { 103 return NaN() 104 } 105 if yi >= 1<<63 { 106 // yi is a large even int that will lead to overflow (or underflow to 0) 107 // for all x except -1 (x == 1 was handled earlier) 108 switch { 109 case x == -1: 110 return 1 111 case (Abs(x) < 1) == (y > 0): 112 return 0 113 default: 114 return Inf(1) 115 } 116 } 117 118 // ans = a1 * 2**ae (= 1 for now). 119 a1 := 1.0 120 ae := 0 121 122 // ans *= x**yf 123 if yf != 0 { 124 if yf > 0.5 { 125 yf-- 126 yi++ 127 } 128 a1 = Exp(yf * Log(x)) 129 } 130 131 // ans *= x**yi 132 // by multiplying in successive squarings 133 // of x according to bits of yi. 134 // accumulate powers of two into exp. 135 x1, xe := Frexp(x) 136 for i := int64(yi); i != 0; i >>= 1 { 137 if xe < -1<<12 || 1<<12 < xe { 138 // catch xe before it overflows the left shift below 139 // Since i !=0 it has at least one bit still set, so ae will accumulate xe 140 // on at least one more iteration, ae += xe is a lower bound on ae 141 // the lower bound on ae exceeds the size of a float64 exp 142 // so the final call to Ldexp will produce under/overflow (0/Inf) 143 ae += xe 144 break 145 } 146 if i&1 == 1 { 147 a1 *= x1 148 ae += xe 149 } 150 x1 *= x1 151 xe <<= 1 152 if x1 < .5 { 153 x1 += x1 154 xe-- 155 } 156 } 157 158 // ans = a1*2**ae 159 // if y < 0 { ans = 1 / ans } 160 // but in the opposite order 161 if y < 0 { 162 a1 = 1 / a1 163 ae = -ae 164 } 165 return Ldexp(a1, ae) 166} 167