Source file src/pkg/math/pow.go

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

View as plain text