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

[J-F Garnier]

(11-30-2022 07:01 PM)Valentin Albillo Wrote:  

J-F Garnier Wrote:Clearly brute force is not the way to go, but I can't figure out a short-cut.

Stick to it, please, I'm sure you'll eventually manage.

The missing point for me was that there are efficient approximations of the Harmonic numbers, and that we can easily find the point when the partial sum can be reduced to log(2).
Now I think I have an idea to calculate the sum in a different (but maybe equivalent) way than Werner, and my understanding is that it will also prove that the sum is convergent.

Edit: Result 2.08637766501 confirmed I will post my HP-71B program and comments soon...

J-F

[Werner]

I found another, similar way, but I fear it's the same as J-F's, so I will hold out till he posted his solution ;-)
Werner

[Valentin Albillo]

  
Hi, all,

A last round of comments before I post my original solution next Sunday night:

John Keith Wrote:Unfortunately my program was just a brute force one and as many have found, brute force is a dead end for this challenge.

Indeed.

John Keith Wrote:After summing 2^17 terms and seeing no convergence, it was obvious that I was getting nowhere.

It would need summing much more terms than that to get just one roughly correct digit.

Quote:I was never anywhere near an analytical solution such as Werner's. Also as you mentioned, the recursion did not turn out to be a problem. It was an interesting and enjoyable challenge, though. Maybe next time...

Glad you liked it. The analytics aren't that difficult at all, it's just getting the correct idea. Thank you very much for your interest and for doing your best to solve it. Maybe next time, though Problems 4, 5 and 6 are allegedly more difficult (though still solvable using an HP-71B) ... or not !

FLISZT Wrote:About the code panels, I will say that it's a pity that you have to scroll as soon as the code is a bit long.

For me, the problem with CODE panels is that when I create a PDF with the thread (which I always do to upload it to my site), the code within them gets truncated and so becomes useless. Thanks for your participation in this thread, Bruno.

Werner Wrote:Apologies for jumping the gun! But I did wait for a full day..

Most commendable. Well, someone had to be first, it might as well be you.

Werner Wrote:Yes, well. The trick I used only works if the sum is convergent to begin with, something I haven't been able to prove.

You can trust me, it is a convergent series and it's not that hard to prove, surely Albert (Chan ?) can provide you with one proof or seven, else I'll obligue.

Werner Wrote:Using the original asymptotic series [...] (thanks^2, Albert!) [...] now, I need the definition only for n<32 instead of 128, and the running time on a real 42S went down to 36 seconds.

Excellent run time indeed, you should've posted "thanks^3" to Albert instead, and when he provides you with the convergence proof you can up it to "thanks^4"  

J-F Garnier Wrote:Now I think I have an idea to calculate the sum in a different (but maybe equivalent) way than Werner, and my understanding is that it will also prove that the sum is convergent.

Good. See ? Perseveration was the key and eventually you did manage, as I said you would.

J-F Garnier Wrote:Edit: Result 2.08637766501 confirmed. I will post my HP-71B program and comments soon...

Perfect. I can confirm that the result is indeed correct to 12 digits. Eagerly waiting for you to post your HP-71B program and comments. You have till next Sunday night.

Werner Wrote:I found another, similar way, but I fear it's the same as J-F's, so I will hold out till he posted his solution ;-)

Most considerate of you, J-F will surely be most pleased.

V.

[J-F Garnier]

Here is my "solution", based on Werner's reasoning, so most credit due to him:-)

I restarted from the step below, changing the K+1 limit (the point from where we use ln(2) as an approximation of the partial sum) to K. It's a just a notation:

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

Here we can do to the last part S3 = (1/F(K) + 1/F(K+1) + .... ) what we already did for the first part of the sum,
i.e. writing F(K) = K.F(d) with d=number of binary digits of K.
If we choose K=2^(M-1), we have F(K) = K.F(M)

S3 = (H(2^M-1)-H(2^(M-1)-1)/F(M) + ... + (H(2^(K-1)-1)-H(2^(K-2)-1))/F(K-1)
+ LN(2)*(1/F(K) + 1/F(K+1) + .... )

And here we recognize that we can do the same thing again on the last (1/F(K) + 1/F(K+1) + .... ) part, and again and again.
Also the quantity S2= (H(2^M-1)-H(2^(M-1)-1)/F(M) + ... + (H(2^(K-1)-1)-H(2^(K-2)-1))/F(K-1)
has already been computed as part of the beginning of the sum.

So we end with:
S = 1 + 1/2 + 1/6 + (H(7)-H(3))/F(3) + .. + (H(2^(M-1)-1)-H(2^(M-2)-1))/F(M-1) + S2
+ LN(2) * ( S2 + LN(2) * (S2 + LN(2) * (S2 ...

The limit Sx of the quantity ( S2 + LN(2) * ( S2 + LN(2) * (S2 + LN(2) * S2 ... ) is finite and is such as S2 + LN(2) * Sx = Sx
so Sx = S2 / (1 - LN(2))

To compute the whole sum S:
compute S1 = 1 + 1/2 + 1/6 + (H(7)-H(3))/F(3) + .. + (H(2^(M-1)-1)-H(2^(M-2)-1))/F(M-1)
compute S2 = (H(2^M-1)-H(2^(M-1)-1)/F(M) + ... + (H(2^(K-1)-1)-H(2^(K-2)-1))/F(K-1)
compute S = S1 + S2 / (1 - LN(2))

M is chosen such as 2^K=2^(2^(M-1)) >> 1E12 to ensure the partial sum H(2^(K-1)-1)-H(2^(K-2)-1) is accurately represented by ln(2) :
I used M=7 so K=2^6=64.



10 ! SRC12C3
20  ! HP-71 / HP-75 version
30  L=63
40  DIM F(63)
50  !
60  L2=LOG(2) ! used several times
70  ! calculate the F(I)
80  F(1)=1 @ F(2)=2
90  FOR I=3 TO L @ D=INT(LOG(I)/L2)+1 @ F(I)=I*F(D) @ NEXT I
100 !
110 ! compute the sum S0 for K=1..2
120 S0=1/F(1)+1/F(2)+1/F(3)
130 ! compute the sum S1 for K=3..6
140 S1=0
150 FOR K=3 TO 6
160   J1=2^(K-1) @ J2=2*J1-1
170   X=0 @ FOR J=J1 TO J2 @ X=X+1/J @ NEXT J
180   S1=S1+X/F(K)
190 NEXT K
210 ! now compute the sum S2 for K=7..40 using the H approx
230 S2=0
240 FOR K=K TO 40
250   N=2^(K+1) @ N2=N*N
260   X=((((-272/N2+16)/N2-2)/N2+1)/N+1)/N+L2
270   S2=S2+X/F(K)
280 NEXT K
290 ! and complete the sum S2 up to K=L using the LOG(2) approx
310 FOR K=K TO L @ S2=S2+L2/F(K) @ NEXT K
320 !
330 ! now compute the final sum S
340 S=S2/(1-LOG(2))+S1+S0
350 DISP S


Result = 2.08637766501
HP-75D: 9.3s
HP-71B: 13.2s

J-F

[Albert Chan]

Sorry about math sessions. I will keep it short.

Let F(n) = n * F(bits of n), except that F(2)=2, F(1)=1
Let G(n) = sum of n-bits integer reciprocal, except that G(1)=5/4

Let sum = index from 1 to K-1, SUM = index K to infinity

S = sum(1/F) + SUM(1/F)       // definition
S = sum(G/F) + SUM(G/F)      // sum(G/F) accelerated convergence, but still not fast enough

S = sum(G/F) + LN2 * SUM(G/LN2/F)

Since G/LN2 ≥ 1, no matter how big K is, we have:

S ≥ sum(G/F) + LN2 * SUM(1/F)
S ≥ sum(G/F) + LN2 * (S - sum(1/F))

S*(1-LN2) ≥ sum((G-LN2)/F)

No need to hard code conditions for index K, or for G converged to LN2.
Just sum RHS until convergence. It will converge, very quickly.

(G-LN2) part shrink at O(1/2^n), which already can converge without F

[PeterP]

Ok, I was able to get it into the 41 with the below code.
As mentioned earlier I had gotten to the equation with ln(2) on the right side but needed Werner's trick (Thank you!) for getting it solved.
Based on HP41 accuracy, I have split the calculation into three segments:
1) Straight Forward calculation of f(x)
2) Using the approximation H(n) ~ ln(n) + 1/(2n) - 1/(12n^2) + 1/(120n^4)
3) Using the approximation that H(n) - H(n-1) ~ Ln(2)

On the hp 41, the approximation 2) is accurate within the precision of the HP41 from n = 2^5 onwards.

On the hp41, the approximation 3) is accurate within the precision of he HP41 from n = 2^31 onwards.

So I am calculating a straight sum from n = 1 till 2^5-1, then switch to the approximation sum 2) from then onwards until 2^31, after which I use ln(2).

I looked at the approximations and I do believe that they prove that the series converges.

I dont know how to do the nice math font here, so my appologies for the below.

S = sum Valentin asked us to calculate.
f(n) = function that valentin gave us
g(x) = H(2^x-1) - H(2^(x-1)-1)
with H(x) being the Harmonic Series up to x
ap(x) = ln(x) + 1/(2x) - 1/(12x^2) + 1/(120x^4)


S = sum(n = 1 to 2^m-1) of f(n)^-1 + sum(x = m+1 to infinity) of f(x)^-1*g(x)
replacing g(x) with the approximation we note that the approximation is always larger than g(x) as we stop with an addition.

S<=sum(n = 1 to 2^m-1) of f(n)^-1 + sum(x = m+1 to infinity) of f(x)^-1* (ap(2^x-1)-ap(2^(x-1)-1)

We then replace the ap() on the right side with ln(2) after a certain cut off point, noting that Ln(2) is also larger than ap()

S<= sum(n = 1 to 2^m-1) of f(n)^-1 + sum(x = m+1 to p-1) of f(x)^-1* (ap(2^x-1)-ap(2^(x-1)-1) + ln(2) * sum(y=p to infinity) of f(y)^-1

with an infinite sum on the right side, this inequality can only be true if S converges.

Or at least this is how my thinking went.

Runtime on the i41cx emulator is a few seconds, result it produces is 2,088075017. Which means that my assumptions about the correct cross-over points are not correct and I might be able to squeeze out a slightly better result by choosing later cross over points but I am not entirely sure how/why, as the calc can't differentiate between ln(2) and the approximation and the approximation and H(2x-1) - H(x-1) anymore at my current cut off points.

However, my flight has landed and I had given up on this when reading it but then had nothing to do on my flight over (europe to US) and made some progress, then got Werners tip, and was able to finish it on the way back. So I feel pretty happy, and thankful.

Here is the code.
Lbl F calculates Valentin's function
Lbl G calculates the approximation of H

Lbl VA12c
CLA
CLX
STO 10
STO 11
STO 12
2
STO 02
31
STO 00
LBL 00
RCL 00
XEQ F
ST+10
DSE00
GTO 00
31.005
STO 00
LBL 01
RCL 00
INT
XEQ G
STO 11
RCL 00
INT
XEQ F
RCL 11
*
ST+ 12
DSE 00
GTO 01
RCL 02
ln
SIGN
LastX
-
1/x
RCL 12
*
RCL 10
+
BEEP
STOP
LBL F
SIGN
STO M
X<>L
LBL 10
3
x>y?
GTO 12
x<>y
ST*M
LN
RCL 02
LN
/
INT
INCX
GTO 10
LBL 12
X<>Y
RCL M
*
1/x
RTN
LBL G
RCL 02
X<>Y
Y^X
DECX
STO N
RCL 02
LASTx
DECX
Y^X
DECX
STO O
LN
LASTx
RCL 02
*
1/x
+
RCLO
X^2
12
*
1/x
-
RCL O
X^2
x^2
120
*
1/x
+
RCL N
LN
LASTx
RCL 02
*
1/x
+
RCL N
X^2
12
*
1/x
-
RCL N
X^2
X^2
120
*
1/x
+
X<>Y
-
RTN

[ijabbott]

(12-03-2022 01:28 AM)Valentin Albillo Wrote:  

FLISZT Wrote:About the code panels, I will say that it's a pity that you have to scroll as soon as the code is a bit long.

For me, the problem with CODE panels is that when I create a PDF with the thread (which I always do to upload it to my site), the code within them gets truncated and so becomes useless. Thanks for your participation in this thread, Bruno.

Just on that point, the "printable" version of threads does expand the CODE panels, and so does the Lite (Archive) view. They would make your PDF files smaller too! It would be nice if the forum had an option to show the entire thread on one page as that would be even easier to convert to PDF.

[Albert Chan]

(12-03-2022 12:15 PM)Albert Chan Wrote:  Let F(n) = n * F(bits of n), except that F(2)=2, F(1)=1
Let G(n) = sum of n-bits integer reciprocal, except that G(1)=5/4

Implementation details, I do not define F(1) or G(1) for simplicity.
Loops sum Z=(G-LN2)/F, from index of 2, until convergence.

(G(1)-LN2)/F(1) = (5/4-LN2)/1 = 1/4 + (1-LN2)

10 DESTROY ALL @ L2=LN(2) @ SETTIME 0
20 DEF FNB(N)=IP(LN(N+.5)/L2)+1 ! BITS OF INTEGER N
30 DEF FNF(N) @ F=2 @ WHILE N>2 @ F=F*N @ N=FNB(N) @ END WHILE @ FNF=F @ END DEF
40 DEF FNG(N) @ G=0 @ N=2^N-1 @ Y=N*(N-1) ! SUM IN PAIRS
50 FOR X=N+N-1 TO N STEP -4 @ G=G+X/Y @ Y=Y-X-X+4 @ NEXT X @ FNG=G @ END DEF
60 DEF FNZ(N) @ IF N<5 THEN Z=FNG(N)-L2 @ GOTO 80
70 Z=.5^(N+1) @ Z2=Z*Z @ Z=(((-272*Z2+16)*Z2-2)*Z2+1)*Z2+Z
80 FNZ=Z/FNF(N) @ END DEF
100 S=1/4 @ I=1 @ REPEAT @ I=I+1 @ P=S @ S=S+FNZ(I) @ UNTIL P=S
110 DISP S/(1-L2)+1,I,TIME

>run
 2.086377665      31      0.1

Emu/DOS WinXP ≈ 200X --> HP71B runtime about 20 seconds.

[Valentin Albillo]

.
Hi, ijabbot,

Quote:Just on that point, the "printable" version of threads does expand the CODE panels, and so does the Lite (Archive) view. They would make your PDF files smaller too!

Yes, I know that the Printable version does expand the CODE panels but last time I checked, the MathJax (or whatever the name is) mathematical expressions do not appear as properly formated textbook expressions but as the MathJax text, which of course defeats the purpose entirely and looks horrible.

In short, when converted to PDF:

- "normal" version gives you nice math expressions but truncated, useless CODE panels.

- "printable" version does expand the code panels but gives you text MathJax instead of properly formated (graphical) math expressions.

If you know some way to have both (i.e. nice math expression and untruncated CODE panels) I'd be very obligued.

Quote:It would be nice if the forum had an option to show the entire thread on one page as that would be even easier to convert to PDF.

Yep, it would be ideal, at least for me. The old forum did just that, the whole thread in a single page, no matter how many posts in it, which was heaven for conversion to PDF.

I've increased the number of posts per page to the maximum, 50, so that if the thread results in 49 replies posted or less then it'll fit in just one page, i.e., one PDF file. Alas, my previous SRC had 86 posts in all so two pages, two PDF.

If Mr. Hicks would increase the limit of posts per page to 100 instead of 50 the situation would improve greatly.

Thanks for trying to help, have a nice weekend.
Regards.
V.
Edit: Oh, and "Lite version" does this:
Example of Lite version with math expressions, images, and CODE panels.
.

[J-F Garnier]

(12-03-2022 05:30 PM)Albert Chan Wrote:  10 DESTROY ALL @ L2=LN(2) @ SETTIME 0
20 DEF FNB(N)=IP(LN(N+.5)/L2)+1 ! BITS OF INTEGER N
30 DEF FNF(N) @ F=2 @ WHILE N>2 @ F=F*N @ N=FNB(N) @ END WHILE @ FNF=F @ END DEF
40 DEF FNG(N) @ G=0 @ N=2^N-1 @ Y=N*(N-1) ! SUM IN PAIRS
50 FOR X=N+N-1 TO N STEP -4 @ G=G+X/Y @ Y=Y-X-X+4 @ NEXT X @ FNG=G @ END DEF
60 DEF FNZ(N) @ IF N<5 THEN Z=FNG(N)-L2 @ GOTO 80
70 Z=.5^(N+1) @ Z2=Z*Z @ Z=(((-272*Z2+16)*Z2-2)*Z2+1)*Z2+Z
80 FNZ=Z/FNF(N) @ END DEF
100 S=1/4 @ I=1 @ REPEAT @ I=I+1 @ P=S @ S=S+FNZ(I) @ UNTIL P=S
110 DISP S/(1-L2)+1,I,TIME

>run
 2.086377665      31      0.1

Emu/DOS WinXP ≈ 200X --> HP71B runtime about 20 seconds.

Runs in 19.9s on a real HP-71B :-)

J-F

[Werner]

Questions!
- Why does Albert’s code run so much slower than J-F’s, while doing 31 loops vs 40?
- Why is a 71B still so much faster than a 42S for similar programs, having even a slower CPU? (My code is down to 31 secs, and is more like Albert’s than like J-F’s, doing 33 loops).
Cheers, Werner

[J-F Garnier]

(12-04-2022 08:49 AM)Werner Wrote:  Questions!
- Why does Albert’s code run so much slower than J-F’s, while doing 31 loops vs 40?
- Why is a 71B still so much faster than a 42S for similar programs, having even a slower CPU? (My code is down to 31 secs, and is more like Albert’s than like J-F’s, doing 33 loops).
Cheers, Werner

Probably a bit OT question, but especially interesting for me!

- I don't fully understand Albert's code, so I'm not sure if his iterative algorithm can be compared to mine/yours.
But there may be at least two reasons that contribute to slower speed:
use of a user function to get f(n) whereas they are precalculated for me,
and FNx user functions overhead.

- for the relative speed of the 42S and 71B, the 42S indeed is not a fast machine due to the System RPL overhead.
The 32S in pure assembly code (as the 71B with similar clock speed, so comparable), is actually faster than the 42S, this is one of the reasons I prefer it to the 42S. The 32SII is unfortunately not as fast as the 32S, for reasons that are really OT here.
From here, relative execution times (the shorter, the better):
- 4:22 HP-32S Keystroke / RPN
- 5:03 HP-32SII Keystroke / RPN / Ver.0
- 12:12 HP-42S Keystroke / RPN / Ver.C

J-F

[Albert Chan]

(12-04-2022 08:49 AM)Werner Wrote:  - Why does Albert’s code run so much slower than J-F’s, while doing 31 loops vs 40?

> 20 DIM F(63) @ F(2)=2 @ F(3)=6 @ B=4
> 30 FOR X=3 TO 6 @ FOR Y=B TO B+B-1 @ F(Y)=Y*F(X) @ NEXT Y @ B=B+B @ NEXT X
> 80 FNZ=Z/F(N) @ END DEF
>RUN
 2.086377665      31      0.05

Above patch removed FNB(N) and FNF(N), and build list of F(N) instead.

F(N=63) is enough, since (G-LN2) shrink at O(1/2^n)
For 12 decimal digits, we need at most N=12/LOG10(2) ≈ 40
Z = (G-LN2)/F, with F growing faster than primes, O(n*ln(n)), it needed even less than that.

Cached F doubled program speed (translated to HP71B runtime of about 10s)

[Albert Chan]

Hi, J-F Garnier

(12-03-2022 08:54 AM)J-F Garnier Wrote:  To compute the whole sum S:
compute S1 = 1 + 1/2 + 1/6 + (H(7)-H(3))/F(3) + .. + (H(2^(M-1)-1)-H(2^(M-2)-1))/F(M-1)
compute S2 = (H(2^M-1)-H(2^(M-1)-1)/F(M) + ... + (H(2^(K-1)-1)-H(2^(K-2)-1))/F(K-1)
compute S = S1 + S2 / (1 - LN(2))

(12-03-2022 12:15 PM)Albert Chan Wrote:  Let sum = index from 1 to K-1, SUM = index K to infinity

S = sum(1/F) + SUM(1/F)       // definition
S = sum(G/F) + SUM(G/F)      // sum(G/F) accelerated convergence, but still not fast enough
...
S*(1-LN2) ≥ sum((G-LN2)/F)

Both are exactly the same, except mine start from 1, yours start from M ≠ 1
If we re-define sum = index from M to K-1, we have:

(S-S1)*(1-LN2) ≥ sum((G-LN2)/F) = S2

[J-F Garnier]

(12-04-2022 10:16 AM)Albert Chan Wrote:  Both are exactly the same, except mine start from 1, yours start from M ≠ 1
If we re-define sum = index from M to K-1, we have:

(S-S1)*(1-LN2) ≥ sum((G-LN2)/F) = S2

(12-03-2022 12:15 PM)Albert Chan Wrote:  S = sum(G/F) + LN2 * SUM(G/LN2/F)

Since G/LN2 ≥ 1, no matter how big K is, we have:

S ≥ sum(G/F) + LN2 * SUM(1/F)
S ≥ sum(G/F) + LN2 * (S - sum(1/F))

S*(1-LN2) ≥ sum((G-LN2)/F)

No need to hard code conditions for index K, or for G converged to LN2.
Just sum RHS until convergence. It will converge, very quickly.

What I don't understand in your reasoning is how you moved from:
S*(1-LN2) ≥ sum((G-LN2)/F)
to
S*(1-LN2) = sum((G-LN2)/F) when the sum has converged.

J-F
(not a mathematician)

[Albert Chan]

(12-04-2022 11:58 AM)J-F Garnier Wrote:  What I don't understand in your reasoning is how you moved from:
S*(1-LN2) ≥ sum((G-LN2)/F)
to
S*(1-LN2) = sum((G-LN2)/F) when the sum has converged.

sum index is from 1 to K-1, but we never specify what K is.
If RHS converge, it meant K can be as big as we wanted (literally, K = infinity)

Another way to look at this:

Last dropped term contributed ≤ 1/2 ULP (definition of "converged")
Because of (G-LN2) shrink rate of O(1/2^n), next dropped term ≤ 1/4 ULP, next ≤ 1/8 ULP, ...

Adding effect of growing F, maximum sum of all missing terms will never contribute 1 ULP.
Numerically, we can turn inequality into equality.

[PeterP]

Interesting improvement to the HP41 result:

At least on the emulator, the accuracy of ln is not good enough to calculate the number of digits in the straight forward fashion of ln(x)/ln(2). For example, for 16, this results in 3.9999999999, rather than 4.

Adding an ulp to it before taking the int fixes the problem, after which my code delivers the expected 2.086377665.

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

[J-F Garnier]

(12-04-2022 09:34 PM)Valentin Albillo Wrote:  I hope you enjoyed it and even learned a little from it.

Thanks for this problem, Valentin!

Yes, I really enjoyed it and leaned a lot about the harmonic series.
I only had the souvenir that the 1/1+1/2+1/3... series is divergent and behaves about as the log function, and even didn't know the Hn notation for the harmonic numbers.

It was really interesting to see different approaches and reasoning, including yours, and the corresponding different programs, all with similar results and performances.

Also, at least for me, the contributions of the other members were very valuable to propose a solution, so it is somehow a collective work.

J-F

[Werner]

.. and a last, ultimate contribution.
improvements:
- new stopping criterion S+H/F = S - stops on the 42S when C=33 instead of 38 ;-)
- choice between H(n) definition and approximation is made in a machine-independent way, so Free42 now calculates the result to (near?) 34-digit precision ;-) and the exact same program can be used in the 42S. Had to rewrite H(n) so that the definition code used 2 counters so that large n could be handled on Free42
- started the summation from C=2, as in Albert's code. Made it ever so slightly shorter - but I added the 1/4 at the end, not at the beginning, gaining 1 digit of accuracy ;-)
- shortened F(n)

on the 42S, max C=33 and S=2.08637766501 in 36 seconds
Free42, C=105 and S=2.086377665005988716089755856734133
and yes, I can shave off a few more bytes here and there..
eg 1 ENTER LN1+X i.o. 1 2 LN
but that would be less readable. These are not 'mini challenges', after all ;-)

00 { 167-Byte Prgm }
01▸LBL "VA3"
02 2
03 STO "C"
04 CLX
05 STO "S"

06▸LBL 10 @ main loop
07 RCL "C"
08 XEQ H
09 1
10 RCL "C"
11 XEQ F
12 ÷
13 RCL+ "S"
14 ENTER
15 X<> "S"
16 X≠Y?
17 ISG "C"
18 X≠Y? @ aff
19 GTO 10

20 4 @ wrap-up
21 1/X
22 RCL+ "S"
23 1
24 2
25 LN
26 -
27 ÷
28 1
29 +
30 RTN

31▸LBL D @ binary digits of an integer
32 CLA
33 BINM
34 ARCL ST X
35 CLX
36 ALENG
37 EXITALL
38 RTN

39▸LBL 04
40 R↓
41 XEQ D
42▸LBL F @ F(n), call with 1 n
43 STO× ST Y
44 3
45 X≤Y?
46 GTO 04
47 R↓
48 R↓
49 RTN

50▸LBL H @ H(2^N-1) - H(2^(N-1)-1) - LN(2), input N
51 2
52 X<>Y
53 Y^X
54 RCL ST X @ if 1/2^n + (2^n)^-10 = 1/2^n then use approximate formula
55 -9
56 Y^X
57 1
58 STO+ ST Y
59 -
60 X=0?
61 GTO 00
62 R↓
63 50
64 %
65 0
66▸LBL 02 @ 1/(n-1) + 1/(n-2) + .. + 1/(n/2)
67 DSE ST Z
68 RCL ST Z
69 1/X
70 +
71 DSE ST Y
72 GTO 02
73 2
74 LN
75 -
76 RTN
77▸LBL 00 @ H(n-1)-H(n/2-1)-ln(2) ~= 1/(2n) + 1/(4n^2) - 1/(8n^4) + 1/(4n^6) - 17/(16n^8)
78 R↓
79 STO+ ST X
80 X^2
81 LASTX
82 1/X
83 272
84 RCL÷ ST Z
85 16
86 -
87 R^
88 ÷
89 2
90 +
91 R^
92 ÷
93 1
94 -
95 R^
96 ÷
97 -
98 END

Cheers, Werner