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