Commit 3c3e68ba authored by Charles L. Dorian's avatar Charles L. Dorian Committed by Russ Cox

math: use ** for exponentiation in comments

R=rsc
CC=golang-dev
https://golang.org/cl/908041
parent 9b1d6332
......@@ -42,7 +42,7 @@ package math
func Acosh(x float64) float64 {
const (
Ln2 = 6.93147180559945286227e-01 // 0x3FE62E42FEFA39EF
Large = 1 << 28 // 2^28
Large = 1 << 28 // 2**28
)
// TODO(rsc): Remove manual inlining of IsNaN
// when compiler does it for us
......@@ -53,9 +53,9 @@ func Acosh(x float64) float64 {
case x == 1:
return 0
case x >= Large:
return Log(x) + Ln2 // x > 2^28
return Log(x) + Ln2 // x > 2**28
case x > 2:
return Log(2*x - 1/(x+Sqrt(x*x-1))) // 2^28 > x > 2
return Log(2*x - 1/(x+Sqrt(x*x-1))) // 2**28 > x > 2
}
t := x - 1
return Log1p(t + Sqrt(2*t+t*t)) // 2 >= x > 1
......
......@@ -28,7 +28,7 @@ package math
// asinh(x) := x if 1+x*x=1,
// := sign(x)*(log(x)+ln2)) for large |x|, else
// := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
// := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
// := sign(x)*log1p(|x| + x**2/(1 + sqrt(1+x**2)))
//
// Asinh(x) calculates the inverse hyperbolic sine of x.
......@@ -40,8 +40,8 @@ package math
func Asinh(x float64) float64 {
const (
Ln2 = 6.93147180559945286227e-01 // 0x3FE62E42FEFA39EF
NearZero = 1.0 / (1 << 28) // 2^-28
Large = 1 << 28 // 2^28
NearZero = 1.0 / (1 << 28) // 2**-28
Large = 1 << 28 // 2**28
)
// TODO(rsc): Remove manual inlining of IsNaN, IsInf
// when compiler does it for us
......@@ -57,13 +57,13 @@ func Asinh(x float64) float64 {
var temp float64
switch {
case x > Large:
temp = Log(x) + Ln2 // |x| > 2^28
temp = Log(x) + Ln2 // |x| > 2**28
case x > 2:
temp = Log(2*x + 1/(Sqrt(x*x+1)+x)) // 2^28 > |x| > 2.0
temp = Log(2*x + 1/(Sqrt(x*x+1)+x)) // 2**28 > |x| > 2.0
case x < NearZero:
temp = x // |x| < 2^-28
temp = x // |x| < 2**-28
default:
temp = Log1p(x + x*x/(1+Sqrt(1+x*x))) // 2.0 > |x| > 2^-28
temp = Log1p(x + x*x/(1+Sqrt(1+x*x))) // 2.0 > |x| > 2**-28
}
if sign {
temp = -temp
......
......@@ -45,7 +45,7 @@ package math
// Atanh(-1) = -Inf
// Atanh(NaN) = NaN
func Atanh(x float64) float64 {
const NearZero = 1.0 / (1 << 28) // 2^-28
const NearZero = 1.0 / (1 << 28) // 2**-28
// TODO(rsc): Remove manual inlining of IsNaN
// when compiler does it for us
// special cases
......
......@@ -27,11 +27,11 @@ const (
// Max is the largest finite value representable by the type.
// Min is the smallest nonzero value representable by the type.
const (
MaxFloat32 = 3.40282346638528859811704183484516925440e+38 /* 2^127 * (2^24 - 1) / 2^23 */
MinFloat32 = 1.401298464324817070923729583289916131280e-45 /* 1 / 2^(127 - 1 + 23) */
MaxFloat32 = 3.40282346638528859811704183484516925440e+38 /* 2**127 * (2**24 - 1) / 2**23 */
MinFloat32 = 1.401298464324817070923729583289916131280e-45 /* 1 / 2**(127 - 1 + 23) */
MaxFloat64 = 1.797693134862315708145274237317043567981e+308 /* 2^1023 * (2^53 - 1) / 2^52 */
MinFloat64 = 4.940656458412465441765687928682213723651e-324 /* 1 / 2^(1023 - 1 + 52) */
MaxFloat64 = 1.797693134862315708145274237317043567981e+308 /* 2**1023 * (2**53 - 1) / 2**52 */
MinFloat64 = 4.940656458412465441765687928682213723651e-324 /* 1 / 2**(1023 - 1 + 52) */
)
// Integer limit values.
......
......@@ -39,7 +39,7 @@ package math
//
// Method:
// 1. For |x| in [0, 0.84375]
// erf(x) = x + x*R(x^2)
// erf(x) = x + x*R(x**2)
// erfc(x) = 1 - erf(x) if x in [-.84375,0.25]
// = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375]
// where R = P/Q where P is an odd poly of degree 8 and
......@@ -49,7 +49,7 @@ package math
//
//
// Remark. The formula is derived by noting
// erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
// erf(x) = (2/sqrt(pi))*(x - x**3/3 + x**5/10 - x**7/42 + ....)
// and that
// 2/sqrt(pi) = 1.128379167095512573896158903121545171688
// is close to one. The interval is chosen because the fix
......@@ -77,7 +77,7 @@ package math
// erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
// erf(x) = 1 - erfc(x)
// where
// R1(z) = degree 7 poly in z, (z=1/x^2)
// R1(z) = degree 7 poly in z, (z=1/x**2)
// S1(z) = degree 8 poly in z
//
// 4. For x in [1/0.35,28]
......@@ -87,7 +87,7 @@ package math
// erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6, else
// erf(x) = sign(x)*(1.0 - tiny)
// where
// R2(z) = degree 6 poly in z, (z=1/x^2)
// R2(z) = degree 6 poly in z, (z=1/x**2)
// S2(z) = degree 7 poly in z
//
// Note1:
......@@ -99,10 +99,10 @@ package math
// Note2:
// Here 4 and 5 make use of the asymptotic series
// exp(-x*x)
// erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
// erfc(x) ~ ---------- * ( 1 + Poly(1/x**2) )
// x*sqrt(pi)
// We use rational approximation to approximate
// g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
// g(s)=f(1/x**2) = log(erfc(x)*x) - x*x + 0.5625
// Here is the error bound for R1/S1 and R2/S2
// |R1/S1 - f(x)| < 2**(-62.57)
// |R2/S2 - f(x)| < 2**(-61.52)
......@@ -189,7 +189,7 @@ const (
func Erf(x float64) float64 {
const (
VeryTiny = 2.848094538889218e-306 // 0x0080000000000000
Small = 1.0 / (1 << 28) // 2^-28
Small = 1.0 / (1 << 28) // 2**-28
)
// special cases
// TODO(rsc): Remove manual inlining of IsNaN, IsInf
......@@ -209,7 +209,7 @@ func Erf(x float64) float64 {
}
if x < 0.84375 { // |x| < 0.84375
var temp float64
if x < Small { // |x| < 2^-28
if x < Small { // |x| < 2**-28
if x < VeryTiny {
temp = 0.125 * (8.0*x + efx8*x) // avoid underflow
} else {
......@@ -266,7 +266,7 @@ func Erf(x float64) float64 {
// Erfc(-Inf) = 2
// Erfc(NaN) = NaN
func Erfc(x float64) float64 {
const Tiny = 1.0 / (1 << 56) // 2^-56
const Tiny = 1.0 / (1 << 56) // 2**-56
// special cases
// TODO(rsc): Remove manual inlining of IsNaN, IsInf
// when compiler does it for us
......@@ -285,7 +285,7 @@ func Erfc(x float64) float64 {
}
if x < 0.84375 { // |x| < 0.84375
var temp float64
if x < Tiny { // |x| < 2^-56
if x < Tiny { // |x| < 2**-56
temp = x
} else {
z := x * x
......
......@@ -59,7 +59,7 @@ package math
//
// 3. Scale back to obtain exp(x):
// From step 1, we have
// exp(x) = 2^k * exp(r)
// exp(x) = 2**k * exp(r)
//
// Special cases:
// exp(INF) is INF, exp(NaN) is NaN;
......@@ -81,7 +81,7 @@ package math
// compiler will convert from decimal to binary accurately enough
// to produce the hexadecimal values shown.
// Exp returns e^x, the base-e exponential of x.
// Exp returns e**x, the base-e exponential of x.
//
// Special cases are:
// Exp(+Inf) = +Inf
......@@ -101,7 +101,7 @@ func Exp(x float64) float64 {
Overflow = 7.09782712893383973096e+02
Underflow = -7.45133219101941108420e+02
NearZero = 1.0 / (1 << 28) // 2^-28
NearZero = 1.0 / (1 << 28) // 2**-28
)
// TODO(rsc): Remove manual inlining of IsNaN, IsInf
......@@ -140,7 +140,7 @@ func Exp(x float64) float64 {
return Ldexp(y, k)
}
// Exp2 returns 2^x, the base-2 exponential of x.
// Exp2 returns 2**x, the base-2 exponential of x.
//
// Special cases are the same as Exp.
func Exp2(x float64) float64 { return Exp(x * Ln2) }
......@@ -34,13 +34,13 @@ package math
// 2. Approximating expm1(r) by a special rational function on
// the interval [0,0.34658]:
// Since
// r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ...
// r*(exp(r)+1)/(exp(r)-1) = 2+ r**2/6 - r**4/360 + ...
// we define R1(r*r) by
// r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r)
// r*(exp(r)+1)/(exp(r)-1) = 2+ r**2/6 * R1(r*r)
// That is,
// R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
// = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r))
// = 1 - r^2/60 + r^4/2520 - r^6/100800 + ...
// = 1 - r**2/60 + r**4/2520 - r**6/100800 + ...
// We use a special Reme algorithm on [0,0.347] to generate
// a polynomial of degree 5 in r*r to approximate R1. The
// maximum error of this polynomial approximation is bounded
......@@ -79,20 +79,20 @@ package math
// = r - E
// 3. Scale back to obtain expm1(x):
// From step 1, we have
// expm1(x) = either 2^k*[expm1(r)+1] - 1
// = or 2^k*[expm1(r) + (1-2^-k)]
// expm1(x) = either 2**k*[expm1(r)+1] - 1
// = or 2**k*[expm1(r) + (1-2**-k)]
// 4. Implementation notes:
// (A). To save one multiplication, we scale the coefficient Qi
// to Qi*2^i, and replace z by (x^2)/2.
// to Qi*2**i, and replace z by (x**2)/2.
// (B). To achieve maximum accuracy, we compute expm1(x) by
// (i) if x < -56*ln2, return -1.0, (raise inexact if x!=inf)
// (ii) if k=0, return r-E
// (iii) if k=-1, return 0.5*(r-E)-0.5
// (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E)
// else return 1.0+2.0*(r-E);
// (v) if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1)
// (vi) if k <= 20, return 2^k((1-2^-k)-(E-r)), else
// (vii) return 2^k(1-((E+2^-k)-r))
// (v) if (k<-2||k>56) return 2**k(1-(E-r)) - 1 (or exp(x)-1)
// (vi) if k <= 20, return 2**k((1-2**-k)-(E-r)), else
// (vii) return 2**k(1-((E+2**-k)-r))
//
// Special cases:
// expm1(INF) is INF, expm1(NaN) is NaN;
......@@ -114,7 +114,7 @@ package math
// to produce the hexadecimal values shown.
//
// Expm1 returns e^x - 1, the base-e exponential of x minus 1.
// Expm1 returns e**x - 1, the base-e exponential of x minus 1.
// It is more accurate than Exp(x) - 1 when x is near zero.
//
// Special cases are:
......@@ -131,7 +131,7 @@ func Expm1(x float64) float64 {
Ln2Hi = 6.93147180369123816490e-01 // 0x3fe62e42fee00000
Ln2Lo = 1.90821492927058770002e-10 // 0x3dea39ef35793c76
InvLn2 = 1.44269504088896338700e+00 // 0x3ff71547652b82fe
Tiny = 1.0 / (1 << 54) // 2^-54 = 0x3c90000000000000
Tiny = 1.0 / (1 << 54) // 2**-54 = 0x3c90000000000000
// scaled coefficients related to expm1
Q1 = -3.33333333333331316428e-02 // 0xBFA11111111110F4
Q2 = 1.58730158725481460165e-03 // 0x3F5A01A019FE5585
......@@ -194,7 +194,7 @@ func Expm1(x float64) float64 {
}
x = hi - lo
c = (hi - x) - lo
} else if absx < Tiny { // when |x| < 2^-54, return x
} else if absx < Tiny { // when |x| < 2**-54, return x
return x
} else {
k = 0
......@@ -223,12 +223,12 @@ func Expm1(x float64) float64 {
return y - 1
}
if k < 20 {
t := Float64frombits(0x3ff0000000000000 - (0x20000000000000 >> uint(k))) // t=1-2^-k
t := Float64frombits(0x3ff0000000000000 - (0x20000000000000 >> uint(k))) // t=1-2**-k
y := t - (e - x)
y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
return y
}
t := Float64frombits(uint64((0x3ff - k) << 52)) // 2^-k
t := Float64frombits(uint64((0x3ff - k) << 52)) // 2**-k
y := x - (e + t)
y += 1
y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent
......
......@@ -6,7 +6,7 @@ package math
// Frexp breaks f into a normalized fraction
// and an integral power of two.
// It returns frac and exp satisfying f == frac × 2^exp,
// It returns frac and exp satisfying f == frac × 2**exp,
// with the absolute value of frac in the interval [½, 1).
func Frexp(f float64) (frac float64, exp int) {
// TODO(rsc): Remove manual inlining of IsNaN, IsInf
......
......@@ -25,11 +25,11 @@ package math
// __ieee754_j0(x), __ieee754_y0(x)
// Bessel function of the first and second kinds of order zero.
// Method -- j0(x):
// 1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
// 1. For tiny x, we use j0(x) = 1 - x**2/4 + x**4/64 - ...
// 2. Reduce x to |x| since j0(x)=j0(-x), and
// for x in (0,2)
// j0(x) = 1-z/4+ z^2*R0/S0, where z = x*x;
// (precision: |j0-1+z/4-z^2R0/S0 |<2**-63.67 )
// j0(x) = 1-z/4+ z**2*R0/S0, where z = x*x;
// (precision: |j0-1+z/4-z**2R0/S0 |<2**-63.67 )
// for x in (2,inf)
// j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
// where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
......@@ -50,13 +50,13 @@ package math
// Method -- y0(x):
// 1. For x<2.
// Since
// y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
// y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x**2/4 - ...)
// therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
// We use the following function to approximate y0,
// y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
// y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x**2
// where
// U(z) = u00 + u01*z + ... + u06*z^6
// V(z) = 1 + v01*z + ... + v04*z^4
// U(z) = u00 + u01*z + ... + u06*z**6
// V(z) = 1 + v01*z + ... + v04*z**4
// with absolute approximation error bounded by 2**-72.
// Note: For tiny x, U/V = u0 and j0(x)~1, hence
// y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27)
......@@ -232,11 +232,11 @@ func Y0(x float64) float64 {
}
// The asymptotic expansions of pzero is
// 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x.
// 1 - 9/128 s**2 + 11025/98304 s**4 - ..., where s = 1/x.
// For x >= 2, We approximate pzero by
// pzero(x) = 1 + (R/S)
// where R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10
// S = 1 + pS0*s^2 + ... + pS4*s^10
// where R = pR0 + pR1*s**2 + pR2*s**4 + ... + pR5*s**10
// S = 1 + pS0*s**2 + ... + pS4*s**10
// and
// | pzero(x)-1-R/S | <= 2 ** ( -60.26)
......@@ -331,13 +331,13 @@ func pzero(x float64) float64 {
}
// For x >= 8, the asymptotic expansions of qzero is
// -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
// -1/8 s + 75/1024 s**3 - ..., where s = 1/x.
// We approximate pzero by
// qzero(x) = s*(-1.25 + (R/S))
// where R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10
// S = 1 + qS0*s^2 + ... + qS5*s^12
// qzero(x) = s*(-1.25 + (R/S))
// where R = qR0 + qR1*s**2 + qR2*s**4 + ... + qR5*s**10
// S = 1 + qS0*s**2 + ... + qS5*s**12
// and
// | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22)
// | qzero(x)/s +1.25-R/S | <= 2**(-61.22)
// for x in [inf, 8]=1/[0,0.125]
var q0R8 = [6]float64{
......
......@@ -25,7 +25,7 @@ package math
// __ieee754_j1(x), __ieee754_y1(x)
// Bessel function of the first and second kinds of order one.
// Method -- j1(x):
// 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
// 1. For tiny x, we use j1(x) = x/2 - x**3/16 + x**5/384 - ...
// 2. Reduce x to |x| since j1(x)=-j1(-x), and
// for x in (0,2)
// j1(x) = x/2 + x*z*R0/S0, where z = x*x;
......@@ -52,13 +52,13 @@ package math
// 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
// 2. For x<2.
// Since
// y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
// y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x**3-...)
// therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
// We use the following function to approximate y1,
// y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
// y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x**2
// where for x in [0,2] (abs err less than 2**-65.89)
// U(z) = U0[0] + U0[1]*z + ... + U0[4]*z^4
// V(z) = 1 + v0[0]*z + ... + v0[4]*z^5
// U(z) = U0[0] + U0[1]*z + ... + U0[4]*z**4
// V(z) = 1 + v0[0]*z + ... + v0[4]*z**5
// Note: For tiny x, 1/x dominate y1 and hence
// y1(tiny) = -2/pi/tiny, (choose tiny<2**-54)
// 3. For x>=2.
......@@ -225,11 +225,11 @@ func Y1(x float64) float64 {
}
// For x >= 8, the asymptotic expansions of pone is
// 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x.
// 1 + 15/128 s**2 - 4725/2**15 s**4 - ..., where s = 1/x.
// We approximate pone by
// pone(x) = 1 + (R/S)
// where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10
// S = 1 + ps0*s^2 + ... + ps4*s^10
// where R = pr0 + pr1*s**2 + pr2*s**4 + ... + pr5*s**10
// S = 1 + ps0*s**2 + ... + ps4*s**10
// and
// | pone(x)-1-R/S | <= 2**(-60.06)
......@@ -324,11 +324,11 @@ func pone(x float64) float64 {
}
// For x >= 8, the asymptotic expansions of qone is
// 3/8 s - 105/1024 s^3 - ..., where s = 1/x.
// 3/8 s - 105/1024 s**3 - ..., where s = 1/x.
// We approximate qone by
// qone(x) = s*(0.375 + (R/S))
// where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10
// S = 1 + qs1*s^2 + ... + qs6*s^12
// where R = qr1*s**2 + qr2*s**4 + ... + qr5*s**10
// S = 1 + qs1*s**2 + ... + qs6*s**12
// and
// | qone(x)/s -0.375-R/S | <= 2**(-61.13)
......
......@@ -64,7 +64,7 @@ func Jn(n int, x float64) float64 {
case x < -MaxFloat64 || x > MaxFloat64: // IsInf(x, 0):
return 0
}
// J(-n, x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
// J(-n, x) = (-1)**n * J(n, x), J(n, -x) = (-1)**n * J(n, x)
// Thus, J(-n, x) = J(n, -x)
if n == 0 {
......@@ -125,7 +125,7 @@ func Jn(n int, x float64) float64 {
} else {
if x < TwoM29 { // x < 2**-29
// x is tiny, return the first Taylor expansion of J(n,x)
// J(n,x) = 1/n!*(x/2)^n - ...
// J(n,x) = 1/n!*(x/2)**n - ...
if n > 33 { // underflow
b = 0
......@@ -135,13 +135,13 @@ func Jn(n int, x float64) float64 {
a := float64(1)
for i := 2; i <= n; i++ {
a *= float64(i) // a = n!
b *= temp // b = (x/2)^n
b *= temp // b = (x/2)**n
}
b /= a
}
} else {
// use backward recurrence
// x x^2 x^2
// x x**2 x**2
// J(n,x)/J(n-1,x) = ---- ------ ------ .....
// 2n - 2(n+1) - 2(n+2)
//
......@@ -187,7 +187,7 @@ func Jn(n int, x float64) float64 {
}
a := t
b = 1
// estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
// estimate log((2/x)**n*n!) = n*log(2/x)+n*ln(n)
// Hence, if n*(log(2n/x)) > ...
// single 8.8722839355e+01
// double 7.09782712893383973096e+02
......
......@@ -5,7 +5,7 @@
package math
// Ldexp is the inverse of Frexp.
// It returns frac × 2^exp.
// It returns frac × 2**exp.
func Ldexp(frac float64, exp int) float64 {
// TODO(rsc): Remove manual inlining of IsNaN, IsInf
// when compiler does it for us
......
......@@ -39,7 +39,7 @@ package math
// minimum (ymin=1.461632144968362245) to maintain monotonicity.
// On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
// Let z = x-ymin;
// lgamma(x) = -1.214862905358496078218 + z^2*poly(z)
// lgamma(x) = -1.214862905358496078218 + z**2*poly(z)
// poly(z) is a 14 degree polynomial.
// 2. Rational approximation in the primary interval [2,3]
// We use the following approximation:
......
......@@ -27,7 +27,7 @@ package math
//
// Method :
// 1. Argument Reduction: find k and f such that
// x = 2^k * (1+f),
// x = 2**k * (1+f),
// where sqrt(2)/2 < 1+f < sqrt(2) .
//
// 2. Approximation of log(1+f).
......
......@@ -24,7 +24,7 @@ package math
//
// Method :
// 1. Argument Reduction: find k and f such that
// 1+x = 2^k * (1+f),
// 1+x = 2**k * (1+f),
// where sqrt(2)/2 < 1+f < sqrt(2) .
//
// Note. If k=0, then f=x is exact. However, if k!=0, then f
......@@ -96,9 +96,9 @@ func Log1p(x float64) float64 {
const (
Sqrt2M1 = 4.142135623730950488017e-01 // Sqrt(2)-1 = 0x3fda827999fcef34
Sqrt2HalfM1 = -2.928932188134524755992e-01 // Sqrt(2)/2-1 = 0xbfd2bec333018866
Small = 1.0 / (1 << 29) // 2^-29 = 0x3e20000000000000
Tiny = 1.0 / (1 << 54) // 2^-54
Two53 = 1 << 53 // 2^53
Small = 1.0 / (1 << 29) // 2**-29 = 0x3e20000000000000
Tiny = 1.0 / (1 << 54) // 2**-54
Two53 = 1 << 53 // 2**53
Ln2Hi = 6.93147180369123816490e-01 // 3fe62e42fee00000
Ln2Lo = 1.90821492927058770002e-10 // 3dea39ef35793c76
Lp1 = 6.666666666666735130e-01 // 3FE5555555555593
......@@ -131,8 +131,8 @@ func Log1p(x float64) float64 {
var iu uint64
k := 1
if absx < Sqrt2M1 { // |x| < Sqrt(2)-1
if absx < Small { // |x| < 2^-29
if absx < Tiny { // |x| < 2^-54
if absx < Small { // |x| < 2**-29
if absx < Tiny { // |x| < 2**-54
return x
}
return x - x*x*0.5
......@@ -175,7 +175,7 @@ func Log1p(x float64) float64 {
}
hfsq := 0.5 * f * f
var s, R, z float64
if iu == 0 { // |f| < 2^-20
if iu == 0 { // |f| < 2**-20
if f == 0 {
if k == 0 {
return 0
......
......@@ -82,11 +82,11 @@ func Pow(x, y float64) float64 {
return Exp(y * Log(x))
}
// ans = a1 * 2^ae (= 1 for now).
// ans = a1 * 2**ae (= 1 for now).
a1 := float64(1)
ae := 0
// ans *= x^yf
// ans *= x**yf
if yf != 0 {
if yf > 0.5 {
yf--
......@@ -95,7 +95,7 @@ func Pow(x, y float64) float64 {
a1 = Exp(yf * Log(x))
}
// ans *= x^yi
// ans *= x**yi
// by multiplying in successive squarings
// of x according to bits of yi.
// accumulate powers of two into exp.
......@@ -113,7 +113,7 @@ func Pow(x, y float64) float64 {
}
}
// ans = a1*2^ae
// ans = a1*2**ae
// if flip { ans = 1 / ans }
// but in the opposite order
if flip {
......
......@@ -31,8 +31,8 @@ package math
// Bit by bit method using integer arithmetic. (Slow, but portable)
// 1. Normalization
// Scale x to y in [1,4) with even powers of 2:
// find an integer k such that 1 <= (y=x*2^(2k)) < 4, then
// sqrt(x) = 2^k * sqrt(y)
// find an integer k such that 1 <= (y=x*2**(2k)) < 4, then
// sqrt(x) = 2**k * sqrt(y)
// 2. Bit by bit computation
// Let q = sqrt(y) truncated to i bit after binary point (q = 1),
// i 0
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment