Cut the Cards

[Albert Chan]

I added a tiny correction to compensate errors for smaller H.

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

Code:

from mpmath import exp, euler, harmonic, psi
f1 = lambda k: k + 6/(1/k+12) - 1
f2 = lambda k: k - 1/(24*k) - .5
f3 = lambda k: k - k/(24*k*k+2.7) - .5    # new formula

def test(lst):
    print 'n\t err(f1)\t err(f2)\t err(f3)'
    for n in lst:
        h = harmonic(n)
        k = exp(h - euler)
        print n,'\t% .10f\t% .10f\t% .10f' % (n-f1(k), n-f2(k), n-f3(k))

>>> test(range(10) + range(10,60,10))

n         err(f1)            err(f2)             err(f3)
0         0.0031607566       0.0127518672       -0.0067666262
1        -0.0003177730       0.0010957186       -0.0001620997
2        -0.0002588484       0.0002719596       -0.0000171164
3        -0.0001712995       0.0001037256       -0.0000035307
4        -0.0001178786       0.0000497922       -0.0000010529
5        -0.0000852046       0.0000275594       -0.0000003957
6        -0.0000641850       0.0000167992       -0.0000001742
7        -0.0000499820       0.0000109785       -0.0000000859
8        -0.0000399758       0.0000075616       -0.0000000463
9        -0.0000326776       0.0000054263       -0.0000000266
10       -0.0000271983       0.0000040243       -0.0000000162
20       -0.0000076840       0.0000005432       -0.0000000006
30       -0.0000035570       0.0000001651       -0.0000000001
40       -0.0000020419       0.0000000705       -0.0000000000
50       -0.0000013228       0.0000000364       -0.0000000000

If needed, we could apply Newton's method to get better estimate.

>>> h = harmonic(10)                # h = 2.9289682539682538
>>> x = f3(exp(h-euler))            # x = 10.000000016205583
>>> x - (harmonic(x) - h) / psi(1, x+1)
mpf('10.0')