The Museum of HP Calculators

HP Forum Archive 18

[ Return to Index | Top of Index ]

Inverse of Harmonic Numbers (hp 33s)
Message #1 Posted by Gerson W. Barbosa on 24 July 2008, 3:12 p.m.

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.

      
Re: Inverse of Harmonic Numbers (hp 33s)
Message #2 Posted by Egan Ford on 24 July 2008, 3:35 p.m.,
in response to message #1 by Gerson W. Barbosa

Interested.

      
Re: Inverse of Harmonic Numbers (hp 33s)
Message #3 Posted by Walter B on 24 July 2008, 3:55 p.m.,
in response to message #1 by Gerson W. Barbosa

Me too.

Cumprimentos,

Walter

            
Re: Inverse of Harmonic Numbers (hp 33s)
Message #4 Posted by Gerson W. Barbosa on 25 July 2008, 9:50 a.m.,
in response to message #3 by Walter B

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.

                  
Re: Inverse of Harmonic Numbers (hp 33s)
Message #5 Posted by Gerson W. Barbosa on 27 July 2008, 1:30 p.m.,
in response to message #4 by Gerson W. Barbosa

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.


[ Return to Index | Top of Index ]

Go back to the main exhibit hall