Inverse of Harmonic Numbers (hp 33s)



#6

N0001 LBL N			M0020 1/x
N0002 11 M0021 +/-
N0003 x<>y M0022 XEQ w
N0004 - M0023 -3
N0005 x>0? M0024 *
N0006 GTO M M0025 SQRTx
N0007 Rv M0026 1/x
N0008 LASTx M0027 1
N0009 5.77215664902E-1 M0028 -
N0010 - M0029 0.5
N0011 e^x M0030 *
N0012 0.5 M0031 RTN
N0013 - W0001 LBL W
N0014 RTN W0002 x<>y
M0001 LBL M W0003 STO A
M0002 Rv W0004 x<>y
M0003 LASTx W0005 STO B
M0004 ENTER W0006 1
M0005 ENTER W0007 +
M0006 4 W0008 LN
M0007 * W0009 STO X
M0008 e^x W0010 FN= F
M0009 150/11 W0011 SOLVE X
M0010 * W0012 RCL A
M0011 1/x W0013 x<>y
M0012 + W0014 RTN
M0013 5.77215664902E-1 F0001 LBL F
M0014 - F0002 RCL X
M0015 2 F0003 e^x
M0016 * F0004 LASTx
M0017 e^x F0005 *
M0018 12 F0006 RCL- B
M0019 * F0007 RTN

PGM LN CK
----------------
N 78 6476
M 189 63A7
W 54 C232
F 21 3529

This will meet Valentin's requirement in
S&SMC#18 (#3) (for N > 3 and N < 1151.8. By the way, 4.012007 returns 30.5235952258 in less than 2 seconds. Formula and remarks later, if someone's interested.

Regards,

Gerson.

Edited: 24 July 2008, 3:14 p.m.


#7

Interested.

#8

Me too.

Cumprimentos,

Walter


#9

Thanks Egan and Walter for your interest!

By analyzing the difference between the actual harmonic number and the well known asymptotic approximation, Hn = ln(n) + EulerGamma, I obtained the following approximation:

                                         1         1          1           1
SUM(k=1,n,1/k) ~ ln(n) + EulerGamma + ----- - -------- + --------- - --------
2 4 6
2 n 12 n 120 n 252 n

Nothing new here though, that's just equation (13) in this mathworld article:

http://mathworld.wolfram.com/HarmonicNumber.html

In order to obtain the inverse function one could think of just isolating n in the expression above. The problem is this cannot be done, at least with usual algebra. The best I was able to do was isolating n in h ~ EulerGamma + ln(n) + 1/(2n) and getting

1
n ~ - --------------------------------

/ 1 \
2 W | - ------------------- |
| h - EulerGamma |
\ 2 e /

where W() is Lambert's W function and EulerGamma is the Euler-Mascheroni constant, 0.577215664901532860606512090082...

Then I added x as a correction term to the exponent of e, isolated it, which was not hard to do, and built up a table of the differences in relation to the actual values and tried some curve fitting. But I did not obtain anything worthwile.

Later I found a better approximation formula in http://groups.google.com/group/sci.math/browse_frm/thread/f2a9c879857e2ac9:

n ~ 1/2*(1/Sqrt(-3*W(-1/(12*e^(2*(h - EulerGamma))))) - 1)
By using the same technique I had tried earlier, I eventually obtained

n ~ 1/2*(1/Sqrt(-3*W(-1/(12*e^(2*(h - EulerGamma + 1/(150/11*e^(4*h))))))) - 1)

The actual constants I obtained were 3.99989985 and 13.63604051 which are very close to 4 and 150/11. This makes me believe the exact constants might have a true mathematical meaning, but I may be wrong.

The approximation works nicely as it can be seen in the table below:

  h                 n             (actual n) 
------------------------------------------------

1.00 1.00022409778 (1)
1.50 2.00002207370 (2)
2.45 6.00000020420 (6)
2.93 10.0108471094 (10.0108470927)
3.50 18.0907442948 (18.0907442943)
4.00 30.1532900556 (30.1532900558)

For h greater than 11, the asymptotic approximation n ~ e^(h - EulerGamma) - 1/2 is used in the program.

The correction term I added will prevent the formula to be evaluated for h greater than 289 or so on the hp 33s due to overflow, but that is not necessary as we have seen. If such were the case, the following interesting property I have observed could be used:

            n(h+k)
lim -------- = e^k
h -> inf n(h)

Regards,

Gerson.

------------------------------------------------------
Harmonic Numbers Program (hp 33s)

A0001 LBL A A0016 252
A0002 LN A0017 /
A0003 LASTx A0018 120
A0004 2 A0019 1/x
A0005 * A0020 -
A0006 1/x A0021 *
A0007 + A0022 12
A0008 STO A A0023 1/x
A0009 LASTx A0024 +
A0010 x^2 A0025 *
A0011 4 A0026 STO- A
A0012 * A0027 RCL A
A0013 ENTER A0028 5.77215664902E-1
A0014 ENTER A0029 +
A0015 ENTER A0030 RTN

LENGTH: 162
CHKSUM: 9920


Edited: 25 July 2008, 12:30 p.m.


#10

Update:

n ~ 1/2*(1/sqrt(-3*W(-1/(12*e^(2*(h-EulerGamma+1/(1649/121*e^(4*h))-1/(1575/113*e^(6*h)))))))-1)

or, for better readability:

/ \
1 | 1 |
n ~ --- | ------------------------------------------------------------------------------------- - 1 |
2 | _______________________________________________________________________________ |
| | - - |
| | | -1 | |
| | -3 W | ------------------------------------------------------------------- | |
| | | / 1 1 \ | |
| | | 2 | h - EulerGamma + ------------ - ----------- + ... (?) | | |
| | | | 1649 4h 1575 6h | | |
| \ | | | ----- e ----- e | | |
| \ | | \ 121 113 / | |
| \| | 12 e | |
\ - - /

The approximation is slightly better:
  h                 n             (actual n) 
------------------------------------------------
1.00 0.99994833094 (1)
1.50 1.99999991839 (2)
2.45 6.00000002735 (6)
2.93 10.0108470957 (10.0108470927)
3.50 18.0907442946 (18.0907442943)
4.00 30.1532900559 (30.1532900558)

I am not sure if I have found the right coefficients in the suggested infinite series, but if such a series exists, the powers of e (4h, 6h, 8h?...) appear to be correct.

HP 33s program:

N0001 LBL N M0025 *
N0002 11 M0026 e^x
N0003 x<>y M0027 12
N0004 - M0028 *
N0005 x>0? M0029 1/x
N0006 GTO M M0030 +/-
N0007 Rv M0031 XEQ W
N0008 LASTx M0032 -3
N0009 5.77215664902E-1 M0033 *
N0010 - M0034 SQRTx
N0011 e^x M0035 1/x
N0012 0.5 M0036 1
N0013 - M0037 -
N0014 RTN M0038 0.5
M0001 LBL M M0039 *
M0002 Rv M0040 RTN
M0003 LASTx W0001 LBL W
M0004 ENTER W0002 x<>y
M0005 ENTER W0003 STO A
M0006 4 W0004 x<>y
M0007 * W0005 STO B
M0008 e^x W0006 1
M0009 1649/121 W0007 +
M0010 * W0008 LN
M0011 1/x W0009 STO X
M0012 x<>y W0010 FN= F
M0012 + W0011 SOLVE X
M0014 LASTx W0012 RCL A
M0015 6 W0013 x<>y
M0016 * W0014 RTN
M0017 e^x F0001 LBL F
M0018 1575/113 F0002 RCL X
M0019 * F0003 e^x
W0020 1/x F0004 LASTx
M0021 - F0005 *
M0022 5.77215664902E-1 F0006 RCL- B
M0023 - F0007 RTN
M0024 2

PGM LN CK
----------------
N 78 6476
M 240 5BB0
W 54 C232
F 21 3529


Edited: 15 Aug 2008, 3:00 p.m.


Possibly Related Threads...
Thread Author Replies Views Last Post
  HP Prime: complex numbers in CAS. Alberto Candel 1 564 12-06-2013, 02:36 PM
Last Post: parisse
  [HP Prime] Plots containing complex numbers bug? Chris Pem10 7 1,072 12-05-2013, 07:40 AM
Last Post: cyrille de Brébisson
  comparing numbers on the WP 34S Kiyoshi Akima 7 955 10-19-2013, 09:28 AM
Last Post: walter b
  HP Prime: Operations with Large Numbers Eddie W. Shore 0 296 10-19-2013, 12:24 AM
Last Post: Eddie W. Shore
  HHC 2013 room numbers David Hayden 2 466 09-20-2013, 05:34 PM
Last Post: sjthomas
  [HP-Prime xcas] operations with complex numbers + BUGs + Request CompSystems 9 1,154 09-08-2013, 10:40 PM
Last Post: CompSystems
  TED Talk: Adam Spencer: Why I fell in love with monster prime numbers Les Bell 3 588 09-05-2013, 12:54 PM
Last Post: Ken Shaw
  Irrationality in numbers....the book Matt Agajanian 4 606 08-30-2013, 04:14 PM
Last Post: Matt Agajanian
  Best HP calculator for crunching numbers rpn fan 51 4,068 08-05-2013, 03:17 PM
Last Post: rpn fan
  HP 50G List of Library numbers already in use? Michael Lopez 2 432 08-04-2013, 10:10 PM
Last Post: Michael Lopez

Forum Jump: