Natural logarithm function (base e)
Method
1. Argument Reduction: find k and f such that
x = 2^k * (1+f),
where sqrt(2)/2 < 1+f < sqrt(2) . 2. Approximation of log(1+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 Remes 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) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
(the values of Lg1 to Lg7 are listed in the program)
and
| 2 14 | -58.45
| Lg1*s +...+Lg7*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
log(1+f) = f - s*(f - R) (if f is not too large)
log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) 3. Finally, log(x) = k*ln2 + log(1+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:
log(x) is NaN with signal if x < 0 (including -INF) ;
log(+INF) is +INF; log(0) is -INF with signal;
log(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.
Definition at line 109 of file log.c.
References Lg1, Lg2, Lg3, Lg4, Lg5, Lg6, Lg7, ln2_hi, ln2_lo, two54, and zero.
Referenced by acosh(), acoth(), acsch(), asech(), asinh(), clog(), cpow(), RealNumber::Log(), log10(), RealNumber::Log2(), and log2p().
111 double hfsq,f,s,z,R,w,t1,t2,dk;
118 if (hx < 0x00100000) {
119 if (((hx&0x7fffffff)|lx)==0)
121 if (hx<0)
return (x-x)/
zero;
126 if (hx >= 0x7ff00000)
return x+x;
129 i = (hx+0x95f64)&0x100000;
133 if((0x000fffff&(2+hx))<3) {
142 R = f*f*(0.5-0.33333333333333333*f);
161 if(k==0)
return f-(hfsq-s*(hfsq+R));
165 if(k==0)
return f-s*(f-R);
#define GET_HIGH_WORD(i, d)
Get the more significant 32 bit int from a double.
static const double two54
#define EXTRACT_WORDS(ix0, ix1, d)
Get two 32 bit ints from a double.
static const double ln2_hi
#define SET_HIGH_WORD(d, v)
Set the more significant 32 bits of a double from an int.
static const double ln2_lo
