Copyright 2009 The Go Authors. All rights reserved. Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.

package math

func ( float64) bool {
	,  := Modf()
	return  == 0 && int64()&1 == 1
}
Special cases taken from FreeBSD's /usr/src/lib/msun/src/e_pow.c updated by IEEE Std. 754-2008 "Section 9.2.1 Special values".
Pow returns x**y, the base-x exponential of y. Special cases are (in order): Pow(x, ±0) = 1 for any x Pow(1, y) = 1 for any y Pow(x, 1) = x for any x Pow(NaN, y) = NaN Pow(x, NaN) = NaN Pow(±0, y) = ±Inf for y an odd integer < 0 Pow(±0, -Inf) = +Inf Pow(±0, +Inf) = +0 Pow(±0, y) = +Inf for finite y < 0 and not an odd integer Pow(±0, y) = ±0 for y an odd integer > 0 Pow(±0, y) = +0 for finite y > 0 and not an odd integer Pow(-1, ±Inf) = 1 Pow(x, +Inf) = +Inf for |x| > 1 Pow(x, -Inf) = +0 for |x| > 1 Pow(x, +Inf) = +0 for |x| < 1 Pow(x, -Inf) = +Inf for |x| < 1 Pow(+Inf, y) = +Inf for y > 0 Pow(+Inf, y) = +0 for y < 0 Pow(-Inf, y) = Pow(-0, -y) Pow(x, y) = NaN for finite x < 0 and finite non-integer y
func (,  float64) float64

func (,  float64) float64 {
	switch {
	case  == 0 ||  == 1:
		return 1
	case  == 1:
		return 
	case IsNaN() || IsNaN():
		return NaN()
	case  == 0:
		switch {
		case  < 0:
			if isOddInt() {
				return Copysign(Inf(1), )
			}
			return Inf(1)
		case  > 0:
			if isOddInt() {
				return 
			}
			return 0
		}
	case IsInf(, 0):
		switch {
		case  == -1:
			return 1
		case (Abs() < 1) == IsInf(, 1):
			return 0
		default:
			return Inf(1)
		}
	case IsInf(, 0):
		if IsInf(, -1) {
			return Pow(1/, -) // Pow(-0, -y)
		}
		switch {
		case  < 0:
			return 0
		case  > 0:
			return Inf(1)
		}
	case  == 0.5:
		return Sqrt()
	case  == -0.5:
		return 1 / Sqrt()
	}

	,  := Modf(Abs())
	if  != 0 &&  < 0 {
		return NaN()
	}
yi is a large even int that will lead to overflow (or underflow to 0) for all x except -1 (x == 1 was handled earlier)
		switch {
		case  == -1:
			return 1
		case (Abs() < 1) == ( > 0):
			return 0
		default:
			return Inf(1)
		}
	}
ans = a1 * 2**ae (= 1 for now).
	 := 1.0
	 := 0
ans *= x**yf
	if  != 0 {
		if  > 0.5 {
			--
			++
		}
		 = Exp( * Log())
	}
ans *= x**yi by multiplying in successive squarings of x according to bits of yi. accumulate powers of two into exp.
	,  := Frexp()
	for  := int64();  != 0;  >>= 1 {
catch xe before it overflows the left shift below Since i !=0 it has at least one bit still set, so ae will accumulate xe on at least one more iteration, ae += xe is a lower bound on ae the lower bound on ae exceeds the size of a float64 exp so the final call to Ldexp will produce under/overflow (0/Inf)
			 += 
			break
		}
		if &1 == 1 {
			 *= 
			 += 
		}
		 *= 
		 <<= 1
		if  < .5 {
			 += 
			--
		}
	}
ans = a1*2**ae if y < 0 { ans = 1 / ans } but in the opposite order
	if  < 0 {
		 = 1 / 
		 = -
	}
	return Ldexp(, )