diff options
Diffstat (limited to 'src/pkg/math/expm1.go')
| -rw-r--r-- | src/pkg/math/expm1.go | 237 |
1 files changed, 0 insertions, 237 deletions
diff --git a/src/pkg/math/expm1.go b/src/pkg/math/expm1.go deleted file mode 100644 index 8f56e15cc..000000000 --- a/src/pkg/math/expm1.go +++ /dev/null @@ -1,237 +0,0 @@ -// Copyright 2010 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 - -// The original C code, the long comment, and the constants -// below are from FreeBSD's /usr/src/lib/msun/src/s_expm1.c -// and came with this notice. The go code is a simplified -// version of the original C. -// -// ==================================================== -// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. -// -// Developed at SunPro, a Sun Microsystems, Inc. business. -// Permission to use, copy, modify, and distribute this -// software is freely granted, provided that this notice -// is preserved. -// ==================================================== -// -// expm1(x) -// Returns exp(x)-1, the exponential of x minus 1. -// -// Method -// 1. Argument reduction: -// Given x, find r and integer k such that -// -// x = k*ln2 + r, |r| <= 0.5*ln2 ~ 0.34658 -// -// Here a correction term c will be computed to compensate -// the error in r when rounded to a floating-point number. -// -// 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 + ... -// we define R1(r*r) by -// 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 + ... -// 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 -// by 2**-61. In other words, -// R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5 -// where Q1 = -1.6666666666666567384E-2, -// Q2 = 3.9682539681370365873E-4, -// Q3 = -9.9206344733435987357E-6, -// Q4 = 2.5051361420808517002E-7, -// Q5 = -6.2843505682382617102E-9; -// (where z=r*r, and the values of Q1 to Q5 are listed below) -// with error bounded by -// | 5 | -61 -// | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2 -// | | -// -// expm1(r) = exp(r)-1 is then computed by the following -// specific way which minimize the accumulation rounding error: -// 2 3 -// r r [ 3 - (R1 + R1*r/2) ] -// expm1(r) = r + --- + --- * [--------------------] -// 2 2 [ 6 - r*(3 - R1*r/2) ] -// -// To compensate the error in the argument reduction, we use -// expm1(r+c) = expm1(r) + c + expm1(r)*c -// ~ expm1(r) + c + r*c -// Thus c+r*c will be added in as the correction terms for -// expm1(r+c). Now rearrange the term to avoid optimization -// screw up: -// ( 2 2 ) -// ({ ( r [ R1 - (3 - R1*r/2) ] ) } r ) -// expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- ) -// ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 ) -// ( ) -// -// = 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)] -// 4. Implementation notes: -// (A). To save one multiplication, we scale the coefficient Qi -// 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)) -// -// Special cases: -// expm1(INF) is INF, expm1(NaN) is NaN; -// expm1(-INF) is -1, and -// for finite argument, only expm1(0)=0 is exact. -// -// Accuracy: -// according to an error analysis, the error is always less than -// 1 ulp (unit in the last place). -// -// Misc. info. -// For IEEE double -// if x > 7.09782712893383973096e+02 then expm1(x) overflow -// -// Constants: -// The hexadecimal values are the intended ones for the following -// constants. The decimal values may be used, provided that the -// compiler will convert from decimal to binary accurately enough -// to produce the hexadecimal values shown. -// - -// 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: -// Expm1(+Inf) = +Inf -// Expm1(-Inf) = -1 -// Expm1(NaN) = NaN -// Very large values overflow to -1 or +Inf. -func Expm1(x float64) float64 - -func expm1(x float64) float64 { - const ( - Othreshold = 7.09782712893383973096e+02 // 0x40862E42FEFA39EF - Ln2X56 = 3.88162421113569373274e+01 // 0x4043687a9f1af2b1 - Ln2HalfX3 = 1.03972077083991796413e+00 // 0x3ff0a2b23f3bab73 - Ln2Half = 3.46573590279972654709e-01 // 0x3fd62e42fefa39ef - Ln2Hi = 6.93147180369123816490e-01 // 0x3fe62e42fee00000 - Ln2Lo = 1.90821492927058770002e-10 // 0x3dea39ef35793c76 - InvLn2 = 1.44269504088896338700e+00 // 0x3ff71547652b82fe - Tiny = 1.0 / (1 << 54) // 2**-54 = 0x3c90000000000000 - // scaled coefficients related to expm1 - Q1 = -3.33333333333331316428e-02 // 0xBFA11111111110F4 - Q2 = 1.58730158725481460165e-03 // 0x3F5A01A019FE5585 - Q3 = -7.93650757867487942473e-05 // 0xBF14CE199EAADBB7 - Q4 = 4.00821782732936239552e-06 // 0x3ED0CFCA86E65239 - Q5 = -2.01099218183624371326e-07 // 0xBE8AFDB76E09C32D - ) - - // special cases - switch { - case IsInf(x, 1) || IsNaN(x): - return x - case IsInf(x, -1): - return -1 - } - - absx := x - sign := false - if x < 0 { - absx = -absx - sign = true - } - - // filter out huge argument - if absx >= Ln2X56 { // if |x| >= 56 * ln2 - if absx >= Othreshold { // if |x| >= 709.78... - return Inf(1) // overflow - } - if sign { - return -1 // x < -56*ln2, return -1.0 - } - } - - // argument reduction - var c float64 - var k int - if absx > Ln2Half { // if |x| > 0.5 * ln2 - var hi, lo float64 - if absx < Ln2HalfX3 { // and |x| < 1.5 * ln2 - if !sign { - hi = x - Ln2Hi - lo = Ln2Lo - k = 1 - } else { - hi = x + Ln2Hi - lo = -Ln2Lo - k = -1 - } - } else { - if !sign { - k = int(InvLn2*x + 0.5) - } else { - k = int(InvLn2*x - 0.5) - } - t := float64(k) - hi = x - t*Ln2Hi // t * Ln2Hi is exact here - lo = t * Ln2Lo - } - x = hi - lo - c = (hi - x) - lo - } else if absx < Tiny { // when |x| < 2**-54, return x - return x - } else { - k = 0 - } - - // x is now in primary range - hfx := 0.5 * x - hxs := x * hfx - r1 := 1 + hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5)))) - t := 3 - r1*hfx - e := hxs * ((r1 - t) / (6.0 - x*t)) - if k != 0 { - e = (x*(e-c) - c) - e -= hxs - switch { - case k == -1: - return 0.5*(x-e) - 0.5 - case k == 1: - if x < -0.25 { - return -2 * (e - (x + 0.5)) - } - return 1 + 2*(x-e) - case k <= -2 || k > 56: // suffice to return exp(x)-1 - y := 1 - (e - x) - y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent - return y - 1 - } - if k < 20 { - 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 - y := x - (e + t) - y += 1 - y = Float64frombits(Float64bits(y) + uint64(k)<<52) // add k to y's exponent - return y - } - return x - (x*e - hxs) // c is 0 -} |