Small RPN exponential routine

[Thomas Klemm]

(01-27-2015 12:56 PM)Eddie W. Shore Wrote:  Is this the origin for the LPN1 function? (I need to be here more often)

Not sure if you mean LN1+X for the HP-41C or log1p in Python:

Quote: log1p(x)

Return the natural logarithm of 1+x (base e).
The result is computed in a way which is accurate for x near zero.

At least in case of HP-calculators the CORDIC algorithm is used and not the Taylor Series of \(\ln(1+x)\):



Thus I might disappoint you: these are not the same.

Cheers
Thomas

[Albert Chan]

(01-27-2015 08:52 PM)Thomas Klemm Wrote:  At least in case of HP-calculators the CORDIC algorithm is used and not the Taylor Series of \(\ln(1+x)\):

There is a faster series for ln(1+x), let y = x/(x+2), thus |y| < 1

ln(1+x)
= ln( ((x+2)+x) / ((x+2)-x) )
= ln( (1 + y) / (1 - y) )
= 2 * (y + y^3/3 + y^5/5 + y^7/7 + ...)

Using post 4 example of ln(0.52) :
For 12 digits accuracy, original ln(1+x) series required 33 terms, y transformed series only 12.

x = 0.52-1 = -0.48
y = -0.48/1.52 = -0.315789473684
ln(0.52) = -0.653926467406 (sumed 12 terms, upto y^21/21)

But, we can do better! |y| can be reduced furthur by scaling.
With 1+x between 1 and exp(1), worst case is if scaling does not help.
-> 1+x <= 2/(1+exp(-1)) = 1.46212
->|y| <= (exp(1)-1)/(3*exp(1)+1) = 0.187691
-> For 12 digits accuracy, term y^15/15 and beyond can be ignored.

Redo above with scaling, sum maximum of 6 terms:
Scaled argument = [1/1.46212, 1.46212) = [0.68394, 1.46212)

ln(0.52) = ln(0.52 * exp(1)) - 1 = ln(1.41350655080) - 1

x = 0.41350655080
y = x/(x+2) = 0.171330196167
ln(0.52) = 2*y + correction - 1 = 0.342660392334 + 0.00341314025930 - 1 = -0.653926467407