A (quite) accurate ln1+x function, or "how close can you get" part II

04092014, 06:44 PM
(This post was last modified: 04092014 07:28 PM by Dieter.)
Post: #1




A (quite) accurate ln1+x function, or "how close can you get" part II
Over the last days there has been some discussion regarding the accurate evaluation of the TVM equation. For small interest rates a dedicated ln1+x function allows significantly better accuracy than the standard method. But such a function is missing on many calculators, as well as its counterpart, a special e^{x}–1 function.
There have been some suggestions on how to emulate a sufficiently accurate ln1+x function on calculators that do not offer one in their function set. On the one hand there is the widely know classic approach as suggested by W. Kahan in the HP15C Advanced Functions Handbook, on the other hand some other ways based on hyperbolic functions have been suggested, both for ln1+x and e^{x}–1. I did some accuracy tests with these methods, the results have been posted in another thread. Among 100.000 random numbers between 1 and 10^{–15} these methods showed errors of about 5...9 units in the last place. I wanted to know if this can be improved, so I tried a new approach. It does not require any exotic hyperbolics and is based on a Taylor series: Let \(u = 1+x\) rounded Then \(ln(1+x) \simeq ln u  \frac{(u1)  x}{u}\) I did a test with this method, using the WP34s emulator with 16 digit standard precision. The following program was used to generate 100.000 random numbers between 1 and 10^{–16}. Dependig on Flag A, either the classic HP/Kahan method (Flag A set) or the new method (Flag A clear) is used. Code: 001 LBL D And here are the results: Select HP/Kahan method: f [SF] [A] The = symbol appears Start: [D] "Running PrOGrAM"... Result in x and y: largest positive error: +8 ULP largest negative error: –5 ULP This matches the error level reported earlier. Now let's see how the new method compares: Select new method: g [CF] [A] The = symbol disappears Start: [D] "Running PrOGrAM"... Result in x and y: largest positive error: +2 ULP largest negative error: –1 ULP That looks much better. Further tests showed the following error distribution: –2 ULP: 0 –1 ULP: 10553 ±0 ULP: 74440 +1 ULP: 14996 +2 ULP: 11 Edit: Another run with 1 milliion random numbers shows the same pattern: –2 ULP: 0 –1 ULP: 105163 ±0 ULP: 744963 +1 ULP: 149766 +2 ULP: 108 So nearly 99,99% of the results are within ±1 ULP. What do you think? Dieter 

04102014, 04:47 AM
Post: #2




RE: A (quite) accurate ln1+x function, or "how close can you get" part II
I'm curious about monotonic and inverse behaviors.
For a and a+ULP, is (f(a+ULP)  f(a)) positive, zero, or negative. Hopefully they'd all be positive, but it can happen that there's a place where you get a string of zeros where f(a) is changing much slower than a. You shouldn't ever find a negative. The value of the error of ainverse(function(a)) and afunction(inverse(a)) would ideally be always zero, of course, but unless both the function and the inverse are absolutely monotonic that won't happen. 

04112014, 07:01 PM
Post: #3




RE: A (quite) accurate ln1+x function, or "how close can you get" part II
Thank you very much for your reply.
(04102014 04:47 AM)htom trites Wrote: For a and a+ULP, is (f(a+ULP)  f(a)) positive, zero, or negative. Hopefully they'd all be positive I assume you mean: "in this case", i.e. for f(x) = ln(1+x). (04102014 04:47 AM)htom trites Wrote: but it can happen that there's a place where you get a string of zeros where f(a) is changing much slower than a. You shouldn't ever find a negative. In a monotonically increasing function, yes. Hm, what about a test with a milliion random numbers? I did one just out of curiosity. There were no negatives. (04102014 04:47 AM)htom trites Wrote: The value of the error of ainverse(function(a)) and afunction(inverse(a)) would ideally be always zero, of course, but unless both the function and the inverse are absolutely monotonic that won't happen. I won't happen either in real life calculators with limited accuracy. ;) Consider for instance sqrt(x) and its inverse x² with, say, 10 digits: 1,414213562 < sqrt(2) < 1,414213563 1,414213562² = 1,999999998 1,414213563² = 2,000000001 So a – inverse(function(a)) is either 2 ULP low or 1 ULP high, although both f(a) and its inverse are strictly monotonic. I did some more tests of the ln1+x approximation suggested above. There is one weak point for negative x between –9,5 · 10^{n} and –10^{n–1}, where n is the working precision (number of significant digits). Here the suggested approximation is typically 5 ULP off, so in this small interval it's not better than the original HP/Kahan method. Otherwise it seems to work fine. Dieter 

01312019, 07:04 PM
(This post was last modified: 02012019 01:30 PM by Albert Chan.)
Post: #4




RE: A (quite) accurate ln1+x function, or "how close can you get" part II
(04112014 07:01 PM)Dieter Wrote: I did some more tests of the ln1+x approximation suggested above. There is one weak point for negative x between –9,5 · 10^{n} and –10^{n–1}, where n is the working precision (number of significant digits). Here the suggested approximation is typically 5 ULP off, so in this small interval it's not better than the original HP/Kahan method. Otherwise it seems to work fine. Excess ULP error is due to correction *lowering* decimal exponent. It does not limited to the edge of working precision. (note: above exponents had the sign wrong) Example, crossing 0.001 boundary: 0.001 = LN(1  0.0009995001666 ...), so try around the edge, say X = 0.00099950016 LN(1+X) = LN( 0.9990004998 ) = 1.000000033e3 (error ~ 0.4 ulp) correction = (X+11X) / (1+X) = +4.004002001e11 (all digits correct) log1p(X) ~ LN(1+X) + correction = 9.9999999930e4 (error = 4 ULP, exponent down 1) Actual error, either absolute (4e13) or relative (4e10) are not affected. 

02012019, 04:26 PM
Post: #5




RE: A (quite) accurate ln1+x function, or "how close can you get" part II
(01312019 07:04 PM)Albert Chan Wrote: Excess ULP error is due to correction *lowering* decimal exponent. To avoid excess ULP error, we like correction same sign as X Y = 1+X, roundedtoward 1.0 log1p(X) ~ LN(Y)  (Y1X)/Y Previous example, log1p(X = 0.00099950016) : Y = roundtoward1 of 1+X = 0.9990004999 (10 digits) log1p(X) ~ LN(Y)  (Y1X)/Y = 9.999999333e4  6.006003001e11 = 9.999999934e4 (all digits correct) 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)