[VA] SRC #012c - Then and Now: Sum
|
12-04-2022, 09:34 PM
(This post was last modified: 12-06-2022 11:21 PM by Valentin Albillo.)
Post: #38
|
|||
|
|||
RE: [VA] SRC #012c - Then and Now: Sum
Hi, all, Well, a full week has elapsed since I posted my OP, which has already passed the ~1,700 views mark (as expected, when the difficulty increases the number of posts and/or views decreases,) and I've got a number of solutions and/or comments, namely by Werner, J-F Garnier, FLISZT, ThomasF, PeterP, Fernando del Rey, John Keith and Albert Chan. Thank you very much to all of you for your interest and valuable contributions. Now I'll provide my detailed original solution to this Problem 3 but first a couple' comments:
and in this case using sheer brute force is utterly hopeless. As will be seen below, naïvely adding up terms of the series goes nowhere; matter of fact, after summing more than 2. 1090 terms we get not even one roughly correct digit, never mind 10-12, so this time the programmer has indeed to use some insight to succeed. 2) I expected everyone to first try adding up terms from the series, but after seeing that doing so led nowhere I also expected you would then try some acceleration methods and/or extrapolation methods but no one actually did (myself included) ! That said, this is my detailed sleuthing process and resulting original solution: My sleuthing process First of all, I realized that d(n) = number of binary digits of n doesn't need any conversions to base 2 and binary-digit counting because we simply have d(n) = IP(LOG2(N)+1), which requires an accurate LOG2 function (like the one provided by the HP-71B Math ROM), lest the IP could be off by one. Naïvely using LN(N)/LN(2) will occasionally fail, even in Free42 Decimal, despite its 34-digit accuracy. Knowing that, I then wrote a bit of sleuthing HP-71B code to generate the first 16 terms of the series and then to add up to 100, 1,000, 10,000 and 100,000 terms:
20 END DEF 30 DESTROY ALL @ FOR N=1 TO 16 @ DISP FNF(N); @ NEXT N @ DISP 40 FOR K=2 TO 5 @ N=10^K @ S=0 @ FOR I=1 TO N @ S=S+1/FNF(I) @ NEXT I @ DISP N;S @ NEXT K >RUN 1 2 6 24 30 36 42 192 216 240 264 288 312 336 360 480 100 1.87752 1000 1.89281 10000 1.90075 100000 1.90642 Now, computing a few d(n) values, the series looks like this (Sum 1): and we notice that the value of d(n) stays the same between powers of 2, as for n = 4 to 7 we have d(n) = 3, then for n = 8 to 15 we have d(n) = 4, etc. Also, the denominators have factors 4, 5, 6, 7, ..., which reminds me of the harmonic numbers: and using them the series can now be expressed like this (Sum 2): so we need a quick way to compute H(k), which for small values of k can be done by just summing k terms of its definition, and for large values of k can be done using its asymptotic series expansion: which has an error less than 1/(252 k6), so using it for values of k exceeding IROUND((1/252E-12)1/6) = 40, will be enough to get 12-digit accuracy. Let's check it:
To check if it suffices, I quickly wrote this other raw concoction:
20 FOR N=3 TO M @ S=S+(FNH(2^N-1)-FNH(2^(N-1)-1))/FNF(N)) @ NEXT N 30 DISP USING "3D,2X,D.4D";M;S @ NEXT M 40 DEF FNH(N) @ IF N>40 THEN 60 50 T=0 @ FOR J=1 TO N @ T=T+1/J @ NEXT J @ FNH=T @ END 60 FNH=LN(N)+.577215664902+1/(2*N)-1/(12*N*N)+1/(120*N^4) @ END DEF 70 DEF FNF(N) @ IF N<3 THEN FNF=N ELSE FNF=N*FNF(IP(LOG2(N)+1)) >RUN terms sum ---------------- 100 1.9416 200 1.9472 300 1.9486 Now, if we look carefully at the numerator's expression H(2n − 1) − H(2n − 1 − 1) within the summation, we find that it converges to ln(2) as n tends to infinity ...
20 DISP USING "2D,2(2X,Z.7D)";N;FNH(2^N-1)-FNH(2^(N-1)-1);RES-K @ NEXT N 30 DEF FNH(N) @ IF N>40 THEN FNH=1/120/N^4-1/12/N/N+1/2/N+.577215664902+LN(N) @ END 40 T=0 @ FOR J=1 TO N @ T=T+1/J @ NEXT J @ FNH=T >RUN N H(..)-H(..) Diff-ln(2) -------------------------- 5 0.7090162 0.0158690 10 0.6936357 0.0004885 15 0.6931624 0.0000153 20 0.6931477 0.0000005 25 0.6931472 0.0000000 Since the cause of the slow convergence of Sum 2 is that the numerators tend to a constant (and thus the convergence to zero is provided only by the slowly-increasing denominators f(n),) we can try and force the numerators to tend to zero (instead of ln(2)) by subtracting precisely that very constant ln(2) from each, which is accomplished by subtracting from Sum 2 the product of the original Sum 1 times ln(2), like this: and now all we need to do is to isolate S, which is accomplished by just dividing both sides by 1 - ln(2), obtaining this expression (Sum 3) for the sum of the series: and as we have we just need to sum 2-M-1 = 1E-12 → M = IP(-LOG2(1E-12)-1) = 38 terms to achieve 12-digit accuracy (this value could be hardcoded as a "magic constant" but I've opted for calculating it in the initialization.) My original solution At long last, my original solution is thus this 273-byte 6-liner: (LOG2 is from the Math ROM)
2 FOR N=3 TO M @ S=S+(FNH(2^N-1)-FNH(2^(N-1)-1)-K)/FNF(N) @ NEXT N 3 STD @ DISP "Sum, #terms:";S/(1-K);M 4 DEF FNH(N) @ IF N>40 THEN FNH=1/120/N^4-1/12/N/N+1/2/N+.577215664902+LN(N) @ END 5 T=0 @ FOR J=1 TO N @ T=T+1/J @ NEXT J @ FNH=T 6 DEF FNF(N) @ IF N<3 THEN FNF=N ELSE FNF=N*FNF(IP(LOG2(N)+1)) Line 1 does some initialization, in particular it computes the number of terms to sum. Line 2 is a loop which adds up the previously computed number of terms. Line 3 simply displays the final result and the number of terms used. Lines 4 and 5 are the definition of H(n), the n-th Harmonic number. Line 6 is the recursive definition of f(n). Let's run it: >RUN
For timing, execute instead:
which takes 0.27" on my emulator and 34" on a physical HP-71B. The correct sum is 2.08637766500599+, which rounds to 2.08637766501, so we've got 12 correct digits (save 1 ulp) using just 38 terms of Sum 3. Well, that'll be it for now, I hope you enjoyed it and even learned a little from it. I'll post next Problem 4 after Christmas, so as to avoid having to compete with Xmas for your attention. V. Edit: corrected a very subtle typo. All My Articles & other Materials here: Valentin Albillo's HP Collection |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 4 Guest(s)