Log
Return the (natural) logarithm of x
Method :
Argument Reduction: find k and f such that x = 2^k * (1+f), where sqrt(2)/2 < 1+f < sqrt(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 s3 + 2/5 s5 + ....., = 2s + sR 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) ~ Lg1s +Lg2s +Lg3s +Lg4s +Lg5s +Lg6s +Lg7s (the values of Lg1 to Lg7 are listed in the program) and | 2 14 | -58.45 | Lg1s +...+Lg7s - R(z) | <= 2 | | Note that 2s = f - sf = f - hfsq + shfsq, where hfsq = ff/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)
Finally, log(x) = kln2 + log(1+f). = kln2_hi+(f-(hfsq-(s*(hfsq+R)+kln2_lo))) Here ln2 is split into two floating point number: ln2_hi + ln2_lo, where nln2_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.