YourWishes/jerryscript
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Implement missing Math logarithm functions from ES6 (#3617)

Browse Source 
Math.log2, Math.log10, Math.log1p, Math.expm1

C implementation ported from fdlibm

Part of Issue #3568

JerryScript-DCO-1.0-Signed-off-by: Rafal Walczyna r.walczyna@samsung.com
This commit is contained in:
Rafal Walczyna
2020-04-04 02:00:41 +02:00
committed by GitHub
parent c74256ccba
commit e470b13096
13 changed files with 1124 additions and 16 deletions
Download Patch File Download Diff File Expand all files Collapse all files
jerry-libm/expm1.c
+305
View File
@@ -0,0 +1,305 @@
/* Copyright JS Foundation and other contributors, http://js.foundation
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
* This file is based on work under the following copyright and permission
* notice:
*
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
*
* @(#)s_expm1.c 5.1 93/09/24
*/
#include "jerry-libm-internal.h"
/* 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;
* z = r*r,
* 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.
*/
#define one 1.0
#define huge 1.0e+300
#define tiny 1.0e-300
#define o_threshold 7.09782712893383973096e+02 /* 0x40862E42, 0xFEFA39EF */
#define ln2_hi 6.93147180369123816490e-01 /* 0x3fe62e42, 0xfee00000 */
#define ln2_lo 1.90821492927058770002e-10 /* 0x3dea39ef, 0x35793c76 */
#define invln2 1.44269504088896338700e+00 /* 0x3ff71547, 0x652b82fe */
/* Scaled Q's: Qn_here = 2**n * Qn_above, for R(2*z) where z = hxs = x*x/2: */
#define Q1 -3.33333333333331316428e-02 /* BFA11111 111110F4 */
#define Q2 1.58730158725481460165e-03 /* 3F5A01A0 19FE5585 */
#define Q3 -7.93650757867487942473e-05 /* BF14CE19 9EAADBB7 */
#define Q4 4.00821782732936239552e-06 /* 3ED0CFCA 86E65239 */
#define Q5 -2.01099218183624371326e-07 /* BE8AFDB7 6E09C32D */
double
expm1 (double x)
{
double y, hi, lo, c, e, hxs, hfx, r1;
double_accessor t, twopk;
int k, xsb;
unsigned int hx;
hx = __HI (x);
xsb = hx & 0x80000000; /* sign bit of x */
hx &= 0x7fffffff; /* high word of |x| */
/* filter out huge and non-finite argument */
if (hx >= 0x4043687A)
{
/* if |x|>=56*ln2 */
if (hx >= 0x40862E42)
{
/* if |x|>=709.78... */
if (hx >= 0x7ff00000)
{
unsigned int low;
low = __LO (x);
if (((hx & 0xfffff) | low) != 0)
{
/* NaN */
return x + x;
}
else
{
/* exp(+-inf)-1={inf,-1} */
return (xsb == 0) ? x : -1.0;
}
}
if (x > o_threshold)
{
/* overflow */
return huge * huge;
}
}
if (xsb != 0)
{
/* x < -56*ln2, return -1.0 with inexact */
if (x + tiny < 0.0) /* raise inexact */
{
/* return -1 */
return tiny - one;
}
}
}
/* argument reduction */
if (hx > 0x3fd62e42)
{
/* if |x| > 0.5 ln2 */
if (hx < 0x3FF0A2B2)
{
/* and |x| < 1.5 ln2 */
if (xsb == 0)
{
hi = x - ln2_hi;
lo = ln2_lo;
k = 1;
}
else
{
hi = x + ln2_hi;
lo = -ln2_lo;
k = -1;
}
}
else
{
k = (int) (invln2 * x + ((xsb == 0) ? 0.5 : -0.5));
t.dbl = k;
hi = x - t.dbl * ln2_hi; /* t*ln2_hi is exact here */
lo = t.dbl * ln2_lo;
}
x = hi - lo;
c = (hi - x) - lo;
}
else if (hx < 0x3c900000)
{
/* 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 = one + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5))));
t.dbl = 3.0 - r1 * hfx;
e = hxs * ((r1 - t.dbl) / (6.0 - x * t.dbl));
if (k == 0)
{
/* c is 0 */
return x - (x * e - hxs);
}
else
{
twopk.as_int.hi = 0x3ff00000 + ((unsigned int) k << 20); /* 2^k */
twopk.as_int.lo = 0;
e = (x * (e - c) - c);
e -= hxs;
if (k == -1)
{
return 0.5 * (x - e) - 0.5;
}
if (k == 1)
{
if (x < -0.25)
{
return -2.0 * (e - (x + 0.5));
}
else
{
return one + 2.0 * (x - e);
}
}
if ((k <= -2) || (k > 56))
{
/* suffice to return exp(x)-1 */
y = one - (e - x);
if (k == 1024)
{
y = y * 2.0 * 0x1p1023;
}
else
{
y = y * twopk.dbl;
}
return y - one;
}
t.dbl = one;
if (k < 20)
{
t.as_int.hi = (0x3ff00000 - (0x200000 >> k)); /* t=1-2^-k */
y = t.dbl - (e - x);
y = y * twopk.dbl;
}
else
{
t.as_int.hi = ((0x3ff - k) << 20); /* 2^-k */
y = x - (e + t.dbl);
y += one;
y = y * twopk.dbl;
}
}
return y;
} /* expm1 */
#undef one
#undef huge
#undef tiny
#undef o_threshold
#undef ln2_hi
#undef ln2_lo
#undef invln2
#undef Q1
#undef Q2
#undef Q3
#undef Q4
#undef Q5
jerry-libm/include/math.h
+4
View File
@@ -58,7 +58,11 @@ double atan2 (double, double);
/* Exponential and logarithmic functions. */
double exp (double);
double expm1 (double);
double log (double);
double log1p (double);
double log2 (double);
double log10 (double);
/* Power functions. */
double pow (double, double);
jerry-libm/jerry-libm-internal.h
+4
View File
@@ -89,7 +89,11 @@ double sin (double x);
double tan (double x);
double exp (double x);
double expm1 (double x);
double log (double x);
double log1p (double x);
double log2 (double x);
double log10 (double);
double pow (double x, double y);
double sqrt (double x);
jerry-libm/log10.c
+116
View File
@@ -0,0 +1,116 @@
/* Copyright JS Foundation and other contributors, http://js.foundation
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
* This file is based on work under the following copyright and permission
* notice:
*
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
*
* @(#)e_log10.c 1.3 95/01/18
*/
#include "jerry-libm-internal.h"
/* log10(x)
* Return the base 10 logarithm of x
*
* Method :
* Let log10_2hi = leading 40 bits of log10(2) and
* log10_2lo = log10(2) - log10_2hi,
* ivln10 = 1/log(10) rounded.
* Then
* n = ilogb(x),
* if(n<0) n = n+1;
* x = scalbn(x,-n);
* log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))
*
* Note 1:
* To guarantee log10(10**n)=n, where 10**n is normal, the rounding
* mode must set to Round-to-Nearest.
* Note 2:
* [1/log(10)] rounded to 53 bits has error .198 ulps;
* log10 is monotonic at all binary break points.
*
* Special cases:
* log10(x) is NaN with signal if x < 0;
* log10(+INF) is +INF with no signal; log10(0) is -INF with signal;
* log10(NaN) is that NaN with no signal;
* log10(10**N) = N for N=0,1,...,22.
*
* 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.
*/
#define zero 0.0
#define two54 1.80143985094819840000e+16 /* 0x43500000, 0x00000000 */
#define ivln10 4.34294481903251816668e-01 /* 0x3FDBCB7B, 0x1526E50E */
#define log10_2hi 3.01029995663611771306e-01 /* 0x3FD34413, 0x509F6000 */
#define log10_2lo 3.69423907715893078616e-13 /* 0x3D59FEF3, 0x11F12B36 */
double
log10 (double x)
{
double y, z;
int i, k, hx;
unsigned lx;
double_accessor temp;
hx = __HI (x); /* high word of x */
lx = __LO (x); /* low word of x */
k = 0;
if (hx < 0x00100000)
{
/* x < 2**-1022 */
if (((hx & 0x7fffffff) | lx) == 0)
{
/* log(+-0)=-inf */
return -two54 / zero;
}
if (hx < 0)
{
/* log(-#) = NaN */
return (x - x) / zero;
}
k -= 54;
x *= two54; /* subnormal number, scale up x */
hx = __HI (x); /* high word of x */
}
if (hx >= 0x7ff00000)
{
return x + x;
}
k += (hx >> 20) - 1023;
i = ((unsigned) k & 0x80000000) >> 31;
hx = (hx & 0x000fffff) | ((0x3ff - i) << 20);
y = (double) (k + i);
temp.dbl = x;
temp.as_int.hi = hx;
z = y * log10_2lo + ivln10 * log (temp.dbl);
return z + y * log10_2hi;
} /* log10 */
#undef zero
#undef two54
#undef ivln10
#undef log10_2hi
#undef log10_2lo
jerry-libm/log1p.c
+245
View File
@@ -0,0 +1,245 @@
/* Copyright JS Foundation and other contributors, http://js.foundation
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
* This file is based on work under the following copyright and permission
* notice:
*
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
*
* @(#)s_log1p.c 5.1 93/09/24
*/
#include "jerry-libm-internal.h"
/* log1p(x)
* Method :
* 1. Argument Reduction: find k and f such that
* 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
* may not be representable exactly. In that case, a correction
* term is need. Let u=1+x rounded. Let c = (1+x)-u, then
* log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u),
* and add back the correction term c/u.
* (Note: when x > 2**53, one can simply return log(x))
*
* 2. Approximation of log1p(f).
* Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
* = 2s + 2/3 s**3 + 2/5 s**5 + .....,
* = 2s + s*R
* We use a special Reme algorithm on [0,0.1716] to generate
* a polynomial of degree 14 to approximate R The maximum error
* of this polynomial approximation is bounded by 2**-58.45. In
* other words,
* 2 4 6 8 10 12 14
* R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s +Lp6*s +Lp7*s
* (the values of Lp1 to Lp7 are listed in the program)
* and
* | 2 14 | -58.45
* | Lp1*s +...+Lp7*s - R(z) | <= 2
* | |
* Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
* In order to guarantee error in log below 1ulp, we compute log
* by
* log1p(f) = f - (hfsq - s*(hfsq+R)).
*
* 3. Finally, log1p(x) = k*ln2 + log1p(f).
* = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
* Here ln2 is split into two floating point number:
* ln2_hi + ln2_lo,
* where n*ln2_hi is always exact for |n| < 2000.
*
* Special cases:
* log1p(x) is NaN with signal if x < -1 (including -INF) ;
* log1p(+INF) is +INF; log1p(-1) is -INF with signal;
* log1p(NaN) is that NaN with no signal.
*
* Accuracy:
* according to an error analysis, the error is always less than
* 1 ulp (unit in the last place).
*
* 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.
*
* Note: Assuming log() return accurate answer, the following
* algorithm can be used to compute log1p(x) to within a few ULP:
*
* u = 1+x;
* if(u==1.0) return x ; else
* return log(u)*(x/(u-1.0));
*
* See HP-15C Advanced Functions Handbook, p.193.
*/
#define zero 0.0
#define ln2_hi 6.93147180369123816490e-01 /* 3fe62e42 fee00000 */
#define ln2_lo 1.90821492927058770002e-10 /* 3dea39ef 35793c76 */
#define two54 1.80143985094819840000e+16 /* 43500000 00000000 */
#define Lp1 6.666666666666735130e-01 /* 3FE55555 55555593 */
#define Lp2 3.999999999940941908e-01 /* 3FD99999 9997FA04 */
#define Lp3 2.857142874366239149e-01 /* 3FD24924 94229359 */
#define Lp4 2.222219843214978396e-01 /* 3FCC71C5 1D8E78AF */
#define Lp5 1.818357216161805012e-01 /* 3FC74664 96CB03DE */
#define Lp6 1.531383769920937332e-01 /* 3FC39A09 D078C69F */
#define Lp7 1.479819860511658591e-01 /* 3FC2F112 DF3E5244 */
double
log1p (double x)
{
double hfsq, f, c, s, z, R;
double_accessor u;
int k, hx, hu, ax;
hx = __HI (x);
ax = hx & 0x7fffffff;
c = 0;
k = 1;
if (hx < 0x3FDA827A)
{
/* 1+x < sqrt(2)+ */
if (ax >= 0x3ff00000)
{
/* x <= -1.0 */
if (x == -1.0)
{
/* log1p(-1) = +inf */
return -two54 / zero;
}
else
{
/* log1p(x<-1) = NaN */
return NAN;
}
}
if (ax < 0x3e200000)
{ /* |x| < 2**-29 */
if ((two54 + x > zero) /* raise inexact */
&& (ax < 0x3c900000)) /* |x| < 2**-54 */
{
return x;
}
else
{
return x - x * x * 0.5;
}
}
if ((hx > 0) || hx <= ((int) 0xbfd2bec4))
{
/* sqrt(2)/2- <= 1+x < sqrt(2)+ */
k = 0;
f = x;
hu = 1;
}
}
if (hx >= 0x7ff00000)
{
return x + x;
}
if (k != 0)
{
if (hx < 0x43400000)
{
u.dbl = 1.0 + x;
hu = u.as_int.hi;
k = (hu >> 20) - 1023;
c = (k > 0) ? 1.0 - (u.dbl - x) : x - (u.dbl - 1.0); /* correction term */
c /= u.dbl;
}
else
{
u.dbl = x;
hu = u.as_int.hi;
k = (hu >> 20) - 1023;
c = 0;
}
hu &= 0x000fffff;
/*
* The approximation to sqrt(2) used in thresholds is not
* critical. However, the ones used above must give less
* strict bounds than the one here so that the k==0 case is
* never reached from here, since here we have committed to
* using the correction term but don't use it if k==0.
*/
if (hu < 0x6a09e)
{
/* u ~< sqrt(2) */
u.as_int.hi = hu | 0x3ff00000; /* normalize u */
}
else
{
k += 1;
u.as_int.hi = hu | 0x3fe00000; /* normalize u/2 */
hu = (0x00100000 - hu) >> 2;
}
f = u.dbl - 1.0;
}
hfsq = 0.5 * f * f;
if (hu == 0)
{
/* |f| < 2**-20 */
if (f == zero)
{
if (k == 0)
{
return zero;
}
else
{
c += k * ln2_lo;
return k * ln2_hi + c;
}
}
R = hfsq * (1.0 - 0.66666666666666666 * f);
if (k == 0)
{
return f - R;
}
else
{
return k * ln2_hi - ((R - (k * ln2_lo + c)) - f);
}
}
s = f / (2.0 + f);
z = s * s;
R = z * (Lp1 +
z * (Lp2 + z * (Lp3 + z * (Lp4 + z * (Lp5 + z * (Lp6 + z * Lp7))))));
if (k == 0)
{
return f - (hfsq - s * (hfsq + R));
}
else
{
return k * ln2_hi - ((hfsq - (s * (hfsq + R) + (k * ln2_lo + c))) - f);
}
} /* log1p */
#undef zero
#undef ln2_hi
#undef ln2_lo
#undef two54
#undef Lp1
#undef Lp2
#undef Lp3
#undef Lp4
#undef Lp5
#undef Lp6
#undef Lp7
jerry-libm/log2.c
+160
View File
@@ -0,0 +1,160 @@
/* Copyright JS Foundation and other contributors, http://js.foundation
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
* This file is based on work under the following copyright and permission
* notice:
*
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
*
* @(#)e_log2.c 1.3 95/01/18
*/
#include "jerry-libm-internal.h"
/* log2(x)
* Return the base 2 logarithm of x. See e_log.c and k_log.h for most
* comments.
*
* This reduces x to {k, 1+f} exactly as in e_log.c, then calls the kernel,
* then does the combining and scaling steps
* log2(x) = (f - 0.5*f*f + k_log1p(f)) / ln2 + k
* in not-quite-routine extra precision.
*/
#define zero 0.0
#define two54 1.80143985094819840000e+16 /* 0x43500000, 0x00000000 */
#define ivln2hi 1.44269504072144627571e+00 /* 0x3FF71547, 0x65200000 */
#define ivln2lo 1.67517131648865118353e-10 /* 0x3DE705FC, 0x2EEFA200 */
#define Lg1 6.666666666666735130e-01 /* 0x3FE55555, 0x55555593 */
#define Lg2 3.999999999940941908e-01 /* 0x3FD99999, 0x9997FA04 */
#define Lg3 2.857142874366239149e-01 /* 0x3FD24924, 0x94229359 */
#define Lg4 2.222219843214978396e-01 /* 0x3FCC71C5, 0x1D8E78AF */
#define Lg5 1.818357216161805012e-01 /* 0x3FC74664, 0x96CB03DE */
#define Lg6 1.531383769920937332e-01 /* 0x3FC39A09, 0xD078C69F */
#define Lg7 1.479819860511658591e-01 /* 0x3FC2F112, 0xDF3E5244 */
double
log2 (double x)
{
double f, hfsq, hi, lo, r, val_hi, val_lo, w, y;
int i, k, hx;
unsigned int lx;
double_accessor temp;
hx = __HI (x); /* high word of x */
lx = __LO (x); /* low word of x */
k = 0;
if (hx < 0x00100000)
{ /* x < 2**-1022 */
if (((hx & 0x7fffffff) | lx) == 0)
{
return -two54 / zero; /* log(+-0)=-inf */
}
if (hx < 0)
{
return (x - x) / zero; /* log(-#) = NaN */
}
k -= 54;
x *= two54; /* subnormal number, scale up x */
hx = __HI (x); /* high word of x */
}
if (hx >= 0x7ff00000)
{
return x + x;
}
if (hx == 0x3ff00000 && lx == 0)
{
return zero; /* log(1) = +0 */
}
k += (hx >> 20) - 1023;
hx &= 0x000fffff;
i = (hx + 0x95f64) & 0x100000;
temp.dbl = x;
temp.as_int.hi = hx | (i ^ 0x3ff00000); /* normalize x or x/2 */
k += (i >> 20);
y = (double) k;
f = temp.dbl - 1.0;
hfsq = 0.5 * f * f;
double s, z, R, t1, t2;
s = f / (2.0 + f);
z = s * s;
w = z * z;
t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
R = t2 + t1;
r = s * (hfsq + R);
/*
* f-hfsq must (for args near 1) be evaluated in extra precision
* to avoid a large cancellation when x is near sqrt(2) or 1/sqrt(2).
* This is fairly efficient since f-hfsq only depends on f, so can
* be evaluated in parallel with R. Not combining hfsq with R also
* keeps R small (though not as small as a true `lo' term would be),
* so that extra precision is not needed for terms involving R.
*
* Compiler bugs involving extra precision used to break Dekker's
* theorem for spitting f-hfsq as hi+lo, unless double_t was used
* or the multi-precision calculations were avoided when double_t
* has extra precision. These problems are now automatically
* avoided as a side effect of the optimization of combining the
* Dekker splitting step with the clear-low-bits step.
*
* y must (for args near sqrt(2) and 1/sqrt(2)) be added in extra
* precision to avoid a very large cancellation when x is very near
* these values. Unlike the above cancellations, this problem is
* specific to base 2. It is strange that adding +-1 is so much
* harder than adding +-ln2 or +-log10_2.
*
* This uses Dekker's theorem to normalize y+val_hi, so the
* compiler bugs are back in some configurations, sigh. And I
* don't want to used double_t to avoid them, since that gives a
* pessimization and the support for avoiding the pessimization
* is not yet available.
*
* The multi-precision calculations for the multiplications are
* routine.
*/
hi = f - hfsq;
temp.dbl = hi;
temp.as_int.lo = 0;
lo = (f - hi) - hfsq + r;
val_hi = hi * ivln2hi;
val_lo = (lo + hi) * ivln2lo + lo * ivln2hi;
/* spadd(val_hi, val_lo, y), except for not using double_t: */
w = y + val_hi;
val_lo += (y - w) + val_hi;
val_hi = w;
return val_lo + val_hi;
} /* log2 */
#undef zero
#undef two54
#undef ivln2hi
#undef ivln2lo
#undef Lg1
#undef Lg2
#undef Lg3
#undef Lg4
#undef Lg5
#undef Lg6
#undef Lg7