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

[Werner]

[updated a few error in indexes. Results and code remain the same]
Time for my reasoning and code ;-)

sum of 1/f(x) with f(x) = n*f(d(x)) with d(x) = number of bits to represent x in binary. f(1)=1 and f(2)=2, by definition.

To get the number of bits used by an integer I use:

LBL D
CLA
BINM
ARCL ST X
CLX
ALENG
EXITALL
RTN

which, on a 42S, gets the correct result over the range 1..2^36-1, which is far more than we'll need.

then, calling F with 1 x XEQ F will calculate f(x):

LBL F
3
X>Y?
GTO 00
Rv
STOx ST Y
XEQ D
GTO F
LBL 00
Rv
x
RTN

Let's write out a few values of f(x)

1 1
2 2
3 6
4 24 = 4*6
5 30 = 5*6
6 36 = 6*6
7 42 = 7*6
8 192 = 8*24
9 216 = 9*24
10 240 = 10*24
11 264 = 11*24
..

so the sum becomes

S = 1 + 1/2 + 1/6 + 1/24 + 1/30 + 1/36 + 1/42 + 1/192 + 1/216 + 1/240 + 1/264 + 1/288 + 1/312 + 1/336 + 1/360 + 1/480 + ..

Probably everyone will just have tried to sum up these numbers.. getting nowhere.
The numbers go down very slowly.. for instance, the millionth term is 1/6e8 - which in this case of course does not mean we have reached an accuracy of 1/6e8, as it is not an alternating series.
However, we can combine terms. Take the sequence

1/24 + 1/30 + 1/36 + 1/42 = 1/(4*6) + 1/(5*6) + 1/6*6) + 1/(7*6)

because all numbers 2^n till 2^(n+1)-1 are written with the same number of bits so

S = 1 + 1/2 + 1/6 + (1/4 + 1/5 + 1/6 + 1/7)/6 + (1/8 + 1/9 + .. 1/15)/24 + (1/16 + .. + 1/31)/30 + ..

with H(n) the nth Harmonic number = 1 + 1/2 + 1/3 + .. + 1/n

S = 1 + 1/2 + 1/6 + (H(7)-H(3))/F(3) + (H(15)-H(7))/F(4) + ..

so we can sum up (H(n-1)-H(n/2-1))/F(C) with n=2^C, but we are still stuck with basically the same convergence rate, or lack thereof.
However, for very large n, H(n-1)-H(n/2-1) = ln(2) to working precision. On Free42, this happens for n=2^111, on a real 42S it is n=2^38

So, we sum up the calculated (H(n-1)-H(n/2-1))/F(C), n=2^C, for C=3..K, where K+1 is such that H(n-1)-H(n/2-1) = ln(2), then we add the remainder multiplied by ln(2):

S = 1 + 1/2 + 1/6 + (H(7)-H(3))/F(3) + (H(15)-H(7))/F(4) + .. + (H(2^K-1)-H(2^(K-1)-1))/F(K)
+ LN(2)*(1/F(K+1) + 1/F(K+2) + .... )

if we add in LN(2)*(1/F(1) + 1/(F2) + .. 1/F(K) we get S again on the right side:

S = 5/3 + (H(7)-H(3))/(F(3) + (H(15)-H(7))/F(4) + (H(2^K-1)-H(2^(K-1)-1))/F(K)
- LN(2)*(1/F(1) + 1/F(2) + 1/F(3) + 1/F(4) + ..1/F(K))
+ LN(2)*S

or

S = (5/3 - 1.5*LN(2) + (H(7)-H(3)-LN(2))/F(3) + (H(15)-H(7)-LN(2))/F(4) + .. (H(2^K-1)-H(2^(K-1)-1)-LN(2))/F(K))/(1-LN(2));

To calculate the H(n):
- for small n we take the definition: H(n) = 1/n + 1/(n-1) + 1/(n-2) + .. 1/2 + 1
- larger n we take the approximation H(n) = ln(x + 1/(24*x + 3.7/x)) + gamma, x=n+0.5 and gamma = 0.577.. (thanks Albert!)
Since we only need H(n-1)-H(n/2-1), we don't need the Euler-Mascheroni constant.
For the 42S, n=128 results in a difference of 3e-14.
- H(2^K-1)-H(2^(K-1)-1)-LN(2) calculated with the approximation formula is LN(A)-LN(B)-LN(2) or LN(1 + (A-2B)/2B) as A and 2B grow ever closer
The factor 5/3 - 3/2*LN(2) I calculate as (10 - LN(512))/6 to avoid cancellation

I don't hardcode the value of K, but test when H(2^C-1)-H(2^(C-1)-1) = LN(2)

Running this on Free42 yields 2.0863776650059887.., on the 42S 2.08637766501 in 55 seconds.

00 { 187-Byte Prgm }
01▸LBL "VA3"
02 3
03 STO "C"
04 CLX
05 STO "S"
06▸LBL 10
07 RCL "C"
08 XEQ H
09 X=0?
10 GTO 00
11 1
12 RCL "C"
13 XEQ F
14 ÷
15 STO+ "S"
16 ISG "C"
17 X<>Y
18 GTO 10
19▸LBL 00
20 10
21 512
22 LN
23 -
24 6
25 ÷
26 RCL+ "S"
27 1
28 2
29 LN
30 -
31 ÷
32 RTN
33▸LBL D
34 CLA
35 BINM
36 ARCL ST X
37 CLX
38 ALENG
39 EXITALL
40 RTN
41▸LBL F
42 3
43 X>Y?
44 GTO 00
45 R↓
46 STO× ST Y
47 XEQ D
48 GTO F
49▸LBL 00
50 R↓
51 ×
52 RTN
53▸LBL H @ H(2^N-1) - H(2^(N-1)-1) - LN(2), input N
54 2
55 X<>Y
56 Y^X
57 128 @ 1E3 on Free42
58 X<Y?
59 GTO 00
60 CLX
61 2
62 -
63 0.05
64 %
65 +
66 1
67 + @ n-1,n/2-1
68 0
69▸LBL 02 @ 1/(n-1) + 1/(n-2) + .. + 1/(n/2)
70 RCL ST Y
71 IP
72 1/X
73 +
74 DSE ST Y
75 GTO 02
76 2
77 LN
78 -
79 RTN
80▸LBL 00
81 R↓
82 ENTER
83 XEQ 14
84 X<>Y
85 2
86 ÷
87 XEQ 14
88 STO+ ST X @
89 STO- ST Y
90 ÷
91 LN1+X
92 RTN
93▸LBL 14
94 0.5
95 -
96 3.7
97 RCL÷ ST Y
98 24
99 RCL× ST Z
100 +
101 1/X
102 +
103 END

Cheers, Werner