2008-03-28 14:56:47 -06:00
|
|
|
// 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.
|
|
|
|
|
2008-06-27 18:06:23 -06:00
|
|
|
package math
|
2008-03-28 14:56:47 -06:00
|
|
|
|
|
|
|
|
2009-03-05 14:31:01 -07:00
|
|
|
// Pow returns x**y, the base-x exponential of y.
|
2009-01-20 15:40:40 -07:00
|
|
|
func Pow(x, y float64) float64 {
|
2008-11-20 11:54:02 -07:00
|
|
|
// TODO: x or y NaN, ±Inf, maybe ±0.
|
|
|
|
switch {
|
|
|
|
case y == 0:
|
2009-11-09 13:07:39 -07:00
|
|
|
return 1
|
2008-11-20 11:54:02 -07:00
|
|
|
case y == 1:
|
2009-11-09 13:07:39 -07:00
|
|
|
return x
|
2008-11-20 11:54:02 -07:00
|
|
|
case x == 0 && y > 0:
|
2009-11-09 13:07:39 -07:00
|
|
|
return 0
|
2008-11-20 11:54:02 -07:00
|
|
|
case x == 0 && y < 0:
|
2009-11-09 13:07:39 -07:00
|
|
|
return Inf(1)
|
2008-11-20 11:54:02 -07:00
|
|
|
case y == 0.5:
|
2009-11-09 13:07:39 -07:00
|
|
|
return Sqrt(x)
|
2008-11-20 11:54:02 -07:00
|
|
|
case y == -0.5:
|
2009-11-09 13:07:39 -07:00
|
|
|
return 1/Sqrt(x)
|
2008-11-20 11:54:02 -07:00
|
|
|
}
|
2008-03-28 14:56:47 -06:00
|
|
|
|
2008-11-20 11:54:02 -07:00
|
|
|
absy := y;
|
|
|
|
flip := false;
|
|
|
|
if absy < 0 {
|
|
|
|
absy = -absy;
|
|
|
|
flip = true;
|
|
|
|
}
|
2009-01-22 17:23:44 -07:00
|
|
|
yi, yf := Modf(absy);
|
2008-11-20 11:54:02 -07:00
|
|
|
if yf != 0 && x < 0 {
|
2009-11-09 13:07:39 -07:00
|
|
|
return NaN()
|
2008-11-20 11:54:02 -07:00
|
|
|
}
|
|
|
|
if yi >= 1<<63 {
|
2009-11-09 13:07:39 -07:00
|
|
|
return Exp(y*Log(x))
|
2008-03-28 14:56:47 -06:00
|
|
|
}
|
|
|
|
|
2008-11-20 12:56:48 -07:00
|
|
|
// ans = a1 * 2^ae (= 1 for now).
|
|
|
|
a1 := float64(1);
|
|
|
|
ae := 0;
|
2008-03-28 14:56:47 -06:00
|
|
|
|
2008-11-20 11:54:02 -07:00
|
|
|
// ans *= x^yf
|
|
|
|
if yf != 0 {
|
|
|
|
if yf > 0.5 {
|
|
|
|
yf--;
|
|
|
|
yi++;
|
2008-03-28 14:56:47 -06:00
|
|
|
}
|
2009-10-06 20:40:35 -06:00
|
|
|
a1 = Exp(yf*Log(x));
|
2008-03-28 14:56:47 -06:00
|
|
|
}
|
|
|
|
|
2008-11-20 11:54:02 -07:00
|
|
|
// ans *= x^yi
|
|
|
|
// by multiplying in successive squarings
|
|
|
|
// of x according to bits of yi.
|
|
|
|
// accumulate powers of two into exp.
|
2009-01-22 17:23:44 -07:00
|
|
|
x1, xe := Frexp(x);
|
2008-11-20 12:56:48 -07:00
|
|
|
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--;
|
2008-03-28 14:56:47 -06:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2008-11-20 12:56:48 -07:00
|
|
|
// ans = a1*2^ae
|
2008-11-20 11:54:02 -07:00
|
|
|
// if flip { ans = 1 / ans }
|
|
|
|
// but in the opposite order
|
|
|
|
if flip {
|
2009-10-06 20:40:35 -06:00
|
|
|
a1 = 1/a1;
|
2008-11-20 12:56:48 -07:00
|
|
|
ae = -ae;
|
2008-03-28 14:56:47 -06:00
|
|
|
}
|
2009-01-22 17:23:44 -07:00
|
|
|
return Ldexp(a1, ae);
|
2008-03-28 14:56:47 -06:00
|
|
|
}
|