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
5package math
6
7func 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
38func Pow(x, y float64) float64 {
39	return libc_pow(x, y)
40}
41
42//extern pow
43func libc_pow(float64, float64) float64
44
45func 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