Cut the Cards

[Albert Chan]

(08-26-2020 02:13 PM)Albert Chan Wrote:  I added a tiny correction to compensate errors for smaller H.

HP50g H^-1 :       « 1. Psi + EXP DUPDUP SQ 24. * 2.7 + / - .5 - »

(08-27-2020 10:07 PM)Gerson W. Barbosa Wrote:  How about « 1. Psi + EXP DUP 1. + 4. * INV OVER 1. EXP / + INV OVER 24. * + INV - .5 - » ?

Just to be clear, both corrections are not based on actual continued fraction of e^ψ(x+1/2)
My correction (and probably Gerson's version) are based on curve fitting, for small H.
I picked 2.7 instead of e to get better accuracy at bigger H, it is simply a trade-off.

To build e^ψ(x+1/2) continued fraction with XCas:

XCas> k := exp(Psi(x+0.5))
XCas> limit(k/x, x=inf)        → 1
XCas> limit(1/(k -x)/x, x=inf)        → 24
XCas> limit(1/(1/(k -x) -24x)/x , x=inf)        → 10/37
XCas> limit(1/(1/(1/(k -x) -24x) -10/37*x)/x , x=inf)        → 0 ???

The last one is likely a XCas bug, Mathematica get the coefficient: 689976/74381

\(\qquad\qquad\exp( \psi(x+1/2)) = x
+\frac{1}{24 \cdot x
+\Large\frac{1}{\frac{10}{37} \cdot x
+\frac{1}{\frac{689976}{74381} \cdot x \;
+\; \cdots}}} \)

To confirm, lets try H^-1(H(20)):

XCas> f(x) := 1/(24*x + 1/(10/37*x + 1/(689976/74381*x)))
XCas> h := float(harmonic(20))          → 3.59773965714
XCas> k := exp(h - euler_gamma)      → 20.5020317757
XCas> fsolve(x+f(x) = k, x=k) - 0.5    → 20.0

If f(x) is small, x ≈ k. We can also iterate for x, like this:

XCas> x0 := k
XCas> x0 := k - f(x0)      → 20.5000002012
XCas> x0 := k - f(x0)      → 20.5
XCas> n := x0 - 0.5        → 20.0