The Museum of HP Calculators

HP Forum Archive 18

 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: 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.

Go back to the main exhibit hall