Source file src/math/pow.go

Documentation: 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  	if haveArchPow {
		40  		return archPow(x, y)
		41  	}
		42  	return pow(x, y)
		43  }
		44  
		45  func 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  

View as plain text