Source file src/pkg/math/pow.go

// 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 isOddInt(x float64) bool {
    xi, xf := Modf(x)
    return xf == 0 && int64(xi)&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 Pow(x, y float64) float64 {
    // TODO(rsc): Remove manual inlining of IsNaN, IsInf
    // when compiler does it for us
    switch {
    case y == 0 || x == 1:
        return 1
    case y == 1:
        return x
    case y == 0.5:
        return Sqrt(x)
    case y == -0.5:
        return 1 / Sqrt(x)
    case x != x || y != y: // IsNaN(x) || IsNaN(y):
        return NaN()
    case x == 0:
        switch {
        case y < 0:
            if isOddInt(y) {
                return Copysign(Inf(1), x)
            }
            return Inf(1)
        case y > 0:
            if isOddInt(y) {
                return x
            }
            return 0
        }
    case y > MaxFloat64 || y < -MaxFloat64: // IsInf(y, 0):
        switch {
        case x == -1:
            return 1
        case (Fabs(x) < 1) == IsInf(y, 1):
            return 0
        default:
            return Inf(1)
        }
    case x > MaxFloat64 || x < -MaxFloat64: // IsInf(x, 0):
        if IsInf(x, -1) {
            return Pow(1/x, -y) // Pow(-0, -y)
        }
        switch {
        case y < 0:
            return 0
        case y > 0:
            return Inf(1)
        }
    }

    absy := y
    flip := false
    if absy < 0 {
        absy = -absy
        flip = true
    }
    yi, yf := Modf(absy)
    if yf != 0 && x < 0 {
        return NaN()
    }
    if yi >= 1<<63 {
        return Exp(y * Log(x))
    }

    // ans = a1 * 2**ae (= 1 for now).
    a1 := float64(1)
    ae := 0

    // ans *= x**yf
    if yf != 0 {
        if yf > 0.5 {
            yf--
            yi++
        }
        a1 = Exp(yf * Log(x))
    }

    // ans *= x**yi
    // by multiplying in successive squarings
    // of x according to bits of yi.
    // accumulate powers of two into exp.
    x1, xe := Frexp(x)
    for i := int64(yi); i != 0; i >>= 1 {
        if i&1 == 1 {
            a1 *= x1
            ae += xe
        }
        x1 *= x1
        xe <<= 1
        if x1 < .5 {
            x1 += x1
            xe--
        }
    }

    // ans = a1*2**ae
    // if flip { ans = 1 / ans }
    // but in the opposite order
    if flip {
        a1 = 1 / a1
        ae = -ae
    }
    return Ldexp(a1, ae)
}