[VA] SRC #012c - Then and Now: Sum

[Valentin Albillo]

  
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:

1) As I said in my OP, I quote:


"Some useful advice is to try and find the correct balance between the program doing all the work with no help from you [...] or else using some insight to help speed up the process."

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. 10⁹⁰ 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:

10 DEF FNF(N) @ IF N<3 THEN FNF=N ELSE FNF=N*FNF(IP(LOG2(N)+1))
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

which told me that the series looked like some kind of "weakened" harmonic series, turned convergent but with such a slow convergence rate that finding its value by adding terms was bound to be hopeless, as the sum didn't seem to converge fast enough even when adding up an exponentially increasing number of terms.

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 k⁶), 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:

>FNH(40)+1/41;FNH(41)


4.30293328284    4.30293328284

Now, using n terms of Sum 2 (say, 100 terms) is equivalent to using 2ⁿ - 1 terms of Sum 1 (i.e. 2¹⁰⁰ - 1 = 1.27 . 10³⁰ terms !!,) thus greatly increasing our chances to achieve a decent enough rate of convergence.

To check if it suffices, I quickly wrote this other raw concoction:

10 DESTROY ALL @ FOR M=100 TO 300 STEP 100 @ S=5/3
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

and though 300 terms of this Sum 2 are equivalent to 2³⁰⁰ - 1 = 2.04 . 10⁹⁰ terms of the original Sum 1, it's plainly clear that the convergence rate is still too low, despite the fact that we're adding up exponentially large chunks of the series.

Now, if we look carefully at the numerator's expression H(2ⁿ − 1) − H(2^{n − 1} − 1) within the summation, we find that it converges to ln(2) as n tends to infinity ...

10 DESTROY ALL @ K=LN(2) @ FOR N=5 TO 25 STEP 5
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

... so each numerator tends to the constant ln(2), divided by consecutive values of f(n), and thus after a while it converges no faster than Sum 1, where the constant is 1 instead of ln(2).

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)

1  DESTROY ALL @ K=LN(2) @ M=IP(-LOG2(1E-12)-1) @ S=5/3-3*K/2
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
Sum, #terms: 2.08637766502   38

For timing, execute instead:


>SETTIME 0 @ CALL @ TIME

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.