Szabolcs Nagy
71d23b3103
The log, log2 and log10 functions share a lot of code and to a lesser extent log1p too. A small part of the code was kept separately in __log1p.h, but since it did not capture much of the common code and it was inlined anyway, it did not solve the issue properly. Now the log functions have significant code duplication, which may be resolved later, until then they need to be modified together. logl, log10l, log2l, log1pl: * Fix the sign when the return value should be -inf. * Remove the volatile hack from log10l (seems unnecessary) log1p, log1pf: * Change the handling of small inputs: only |x|<2^-53 is special (then it is enough to return x with the usual subnormal handling) this fixes the sign of log1p(0) in downward rounding. * Do not handle the k==0 case specially (other than skipping the elaborate argument reduction) * Do not handle 1+x close to power-of-two specially (this code was used rarely, did not give much speed up and the precision wasn't better than the general) * Fix the correction term formula (c=1-(u-x) was used incorrectly when x<1 but (double)(x+1)==2, this was not a critical issue) * Use the exact same method for calculating log(1+f) as in log (except in log1p the c correction term is added to the result). log, logf, log10, log10f, log2, log2f: * Use double_t and float_t consistently. * Now the first part of log10 and log2 is identical to log (until the return statement, hopefully this makes maintainence easier). * Most special case formulas were removed (close to power-of-two and k==0 cases), they increase the code size without providing precision or performance benefits (and obfuscate the code). Only x==1 is handled specially so in downward rounding mode the sign of zero is correct (the general formula happens to give -0). * For x==0 instead of -1/0.0 or -two54/0.0, return -1/(x*x) to force raising the exception at runtime. * Arg reduction code is changed (slightly simplified) * The thresholds for arg reduction to [sqrt(2)/2,sqrt(2)] are now consistently the [0x3fe6a09e00000000,0x3ff6a09dffffffff] and the [0x3f3504f3,0x3fb504f2] intervals for double and float reductions respectively (the exact threshold values are not critical) * Remove the obsolete comment for the FLT_EVAL_METHOD!=0 case in log2f (The same code is used for all eval methods now, on i386 slightly simpler code could be used, but we have asm there anyway) all: * Fix signed int arithmetics (using unsigned for bitmanipulation) * Fix various comments
123 lines
3.8 KiB
C
123 lines
3.8 KiB
C
/* origin: FreeBSD /usr/src/lib/msun/src/s_log1p.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.
|
|
* ====================================================
|
|
*/
|
|
/* double log1p(double x)
|
|
* Return the natural logarithm of 1+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 log(1+f): See log.c
|
|
*
|
|
* 3. Finally, log1p(x) = k*ln2 + log(1+f) + c/u. See log.c
|
|
*
|
|
* 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.
|
|
*/
|
|
|
|
#include "libm.h"
|
|
|
|
static const double
|
|
ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
|
|
ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
|
|
Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
|
|
Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
|
|
Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
|
|
Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
|
|
Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
|
|
Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
|
|
Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
|
|
|
|
double log1p(double x)
|
|
{
|
|
union {double f; uint64_t i;} u = {x};
|
|
double_t hfsq,f,c,s,z,R,w,t1,t2,dk;
|
|
uint32_t hx,hu;
|
|
int k;
|
|
|
|
hx = u.i>>32;
|
|
k = 1;
|
|
if (hx < 0x3fda827a || hx>>31) { /* 1+x < sqrt(2)+ */
|
|
if (hx >= 0xbff00000) { /* x <= -1.0 */
|
|
if (x == -1)
|
|
return x/0.0; /* log1p(-1) = -inf */
|
|
return (x-x)/0.0; /* log1p(x<-1) = NaN */
|
|
}
|
|
if (hx<<1 < 0x3ca00000<<1) { /* |x| < 2**-53 */
|
|
/* underflow if subnormal */
|
|
if ((hx&0x7ff00000) == 0)
|
|
FORCE_EVAL((float)x);
|
|
return x;
|
|
}
|
|
if (hx <= 0xbfd2bec4) { /* sqrt(2)/2- <= 1+x < sqrt(2)+ */
|
|
k = 0;
|
|
c = 0;
|
|
f = x;
|
|
}
|
|
} else if (hx >= 0x7ff00000)
|
|
return x;
|
|
if (k) {
|
|
u.f = 1 + x;
|
|
hu = u.i>>32;
|
|
hu += 0x3ff00000 - 0x3fe6a09e;
|
|
k = (int)(hu>>20) - 0x3ff;
|
|
/* correction term ~ log(1+x)-log(u), avoid underflow in c/u */
|
|
if (k < 54) {
|
|
c = k >= 2 ? 1-(u.f-x) : x-(u.f-1);
|
|
c /= u.f;
|
|
} else
|
|
c = 0;
|
|
/* reduce u into [sqrt(2)/2, sqrt(2)] */
|
|
hu = (hu&0x000fffff) + 0x3fe6a09e;
|
|
u.i = (uint64_t)hu<<32 | (u.i&0xffffffff);
|
|
f = u.f - 1;
|
|
}
|
|
hfsq = 0.5*f*f;
|
|
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;
|
|
dk = k;
|
|
return s*(hfsq+R) + (dk*ln2_lo+c) - hfsq + f + dk*ln2_hi;
|
|
}
|