RPN83P: RPN calculator for TI-83+ TI-84+ inspired by HP-42S

[Albert Chan]

(11-22-2023 03:53 PM)Thomas Klemm Wrote:  

(11-21-2023 06:00 PM)bxparks Wrote:   ; This uses ln(1+x) = x * log(1+x) / ((1+x)-1). I think this algorithm is
; faster than above, but apparently it doesn't work on certain computers
; (though I don't know why).
; See https://math.stackexchange.com/questions/175891

What do you mean exactly with "certain computers"?

This was mentioned in Kahan's article,
"How Futile are Mindless Assessments of Roundoff in Floating-Point Computation?", page 22

Kahan Wrote:For example, no spike should mar the graph of
log₂(x) / (x–1) = ( 1 – (x–1)/2 + (x–1)² / 3 – (x–1)³ / 4 + … )/ log(2)
plotted at arguments x near but not 1 . However, here is a spike that has persisted since 1994 in
three M ATLAB versions 4.2 to 6.5 on all my computers (this graph came from an IBM PC)

The quote is for log2, but that is just a scaling factor for log.
My guess is code for log/log2 was bad, causing spikes in plots.

Note that Kahan's formula is actually designed for ratio, log1p(x)/x
(same article (page 19) have formula for expm1(x)/x, with similar idea)

Let y=x+1., z=y-1., error ε = exact - approx = x-z

log1p(x) = log(1+z+ε) = log(1+z) + log(1+ε/(1+z)) ≈ log(1+z) + ε/(1+z)
x = z + ε

If z is small, log(1+z) ≈ z, ε/(1+z) ≈ ε
Numerator and Denominator has about same relative errors, cancelled each other out.

--> log1p(x)/x ≈ log1p(z)/z = log(y)/(y-1)
--> log1p(x) ≈ x * log(y)/(y-1)

I prefer Dieter's formula, with slightly better accuracy: log1p(x) ≈ log(y) - (y-1-x)/y