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

[Valentin Albillo]

  
Hi, all,

After the nice solutions posted for Problem 2 and the 2,500 views mark exceeded (plus no less than three other related threads created by Albert Chan, Thomas Klemm and J-F Garnier,) now's the time for the next part of my SRC #012 - Then and Now, where I'm showing that vintage HP calcs which were able problem-solvers back THEN in the 80's can NOW solve recent tricky problems intended to be tackled with fast modern computers, never mind slow ancient pocket calcs.

In the next weeks I'm proposing six increasingly harder such problems for you to try and solve using your vintage HP calcs while abiding by the mandatory rules summarized here:

You must use VINTAGE HP CALCS (physical/virtual,) coding in either RPN (inc. mcode), RPL (variants existing at the time, inc. SysRPL) or HP-71B languages (inc. BASIC, FORTH, Assembler), so NO XCAS, MATHEMATICA, MAPLE, EXCEL, C/C++/C#, PYTHON, etc., NO LENGTHY MATH SESSIONS and NO CODE PANELS.

On the plus side, you may use any official/popular modules, pacs or libraries available at the time, such as the Math Pac, HP-IL and JPC ROMs for the HP-71B, the Advantage Module, PPC ROM and Extended Memory for the HP-41, and assorted libraries for the RPL models, to name a few.

That said, I'm done with holding back so here you are, a new problem which deals with the sum of an infinite series, namely:

As an aside, you may remember the RPN program featured in my HP Article VA001 - Long Live the HP-11C, which used an efficient algorithm to quickly and accurately compute the sum of alternating series (say 1 - 1/2 + 1/3 - 1/4 + ...), even oscillating (say 1 - 1 + 1 - 1 + ...) or divergent ones (say 1 - 2 + 3 - 4 + ....)

Problem 3:  Sum

Write a sumptuous program to summarily sum this unassuming yet serious series:
    
where f(n) is defined thus:  if n < 3 then f(n) = n else  f(n) = n * f(d(n)), where d(n) = number of binary digits of n.

An example will make it clear; to compute f(10) we proceed as follows:

We need f(10), which is  10 * f(d(10)) = 10 * f(4)   {as 10 = 1010₂ which has 4 binary digits}
     and now we need f(4), which is   4 * f(d(4)) = 4 * f(3)   {as 4 = 100₂ which has 3 binary digits)
          and now we need f(3), which is   3 * f(d(3)) = 3 * f(2)   {as 3 = 11₂ which has 2 binary digits)
               and now we need f(2), which is 2 by definition

and backtracking we have f(3) = 3*2 = 6 and then f(4) = 4*6 = 24 and finally f(10) = 10*24 = 240.

Your program should have no inputs and must output the sum and automatically end. You should strive for 10-12 correct digits and the faster the running time the better.

Some useful advice is to try and find the correct balance between the program doing all the work with no help from you (i.e. sheer brute force, which might take substantial running time,) or else using some insight to help speed up the process. Your choice.

If I see interest, in a week or so I'll post my own original solution for the HP-71B, which is a 6-line program that does the job. In the meantime, let's see your very own clever solutions AND remember the above rules, please.

V.

[FLISZT]

Hi everyone,


This is my very first post on this forum.
I've been lurking around for years, at least 12 if I take into account the old forum.
This "Sum" problem has seemed not too mathematical to me. It inspired me and finally decided me to register.

So, last night I wrote an RPL program on one of my two hp-50g. Then I checked to see if it worked on my trusty hp-28s - released in the 1980s, the 28 series is the only RPL calculator that follows the rules.

To my surprise, the program did not work: I didn't know - or had forgotten - that arithmetic operations on local variables were not yet implemented on these RPL calculators.
After the few necessary corrections, the size of the program increased from 154 to 161 bytes.
The program uses two local variables.

During the day, while rethinking my program, I realized that some instructions were unnecessary and the final size of the program is therefore reduced to 148.5 bytes.

Here are some results to check if my program is OK (?):

f(0) = 0
f(1) = 1
f(2) = 2
f(3) = 6
… which seem obvious at this point!

f(11) = 264
f(12) = 288
f(100 = 4200

… and f(10) = 240, of course!

[Valentin Albillo]

.
Hi, FLISZT,

(11-28-2022 06:19 PM)FLISZT Wrote:  This is my very first post on this forum.
I've been lurking around for years, at least 12 if I take into account the old forum.

Welcome to the forum and congratulations for stopping being a lurker, hope you'll enjoy the place, very nice & knowledgeable people here (mostly).

Quote:This "Sum" problem has seemed not too mathematical to me. It inspired me and finally decided me to register.

Well, thanks for your interest in this thread I recently created and I'm very glad it inspired you to join us, the more we are the merrier. If you like challenges, you might want to have a look at the Challenges section of my HP calcs site (it also features SHARP Pocket Computers).

Quote:So, last night I wrote an RPL program on one of my two hp-50g. Then I checked to see if it worked on my trusty hp-28s [...] To my surprise, the program did not work [...] During the day, while rethinking my program, I realized that some instructions were unnecessary and the final size of the program is therefore reduced to 148.5 bytes.

Good. Please post the code when you feel you're ready to share it with us all.

Quote:Here are some results to check if my program is OK (?):

f(0) = 0, f(1) = 1, f(2) = 2, f(3) = 6 … which seem obvious at this point!

f(11) = 264, f(12) = 288, f(100 = 4200 … and f(10) = 240, of course!

Of course. All of them are correct.

Best regards.
V.

[FLISZT]

Hi Valentin !

Thank you for welcoming me.

Yes I'm ready.

With the 9th Symphony of Bruckner broadcasted at the radio / Eugen Jochum / Bavarian Radio Orchestra (1954) and now the Concerto for Piano & Orch #3 / Robert Casadesus, New York Philh, Dimitri Mitropoulos (live, 1957)… how could it not be the case?!

Here is my attempt:

@ duplication of the value n; m is initialized at 1 (the neutral element of the multiplication…)
<<
DUP 1 -> b m
@ flag 1 is set; BIN mode selected
<<
1 SF
BIN
@ as long as flag 1 is set
WHILE 1 FS?
REPEAT
b 2 > IF
THEN
@ start counting the number of binary digits:
@ conversion of b into a binary object
b R->B
@ conversion of this object into a string (10 ⇒ "# 1010b")
->STR
SIZE
@ size of the string minus 3 characters
@ 2 of them are the "prefix" of the binary object / 1 is the suffix: "b"
3 -
@ n * f(d(n))
DUP
'b' STO
m *
'm' STO
ELSE
@ flag 1 is cleared (signal of ending the "WHILE" loop)
1 CF
END
END
@ The last calculated value of m (or 1… ) is multiplied
@ by the value of n put in the stack at the very beginning
m *
>>
>>


Edit: typo
Edit2 : no code panel!

[Werner]

Well, I have a result..

2.08637766501

I still have to double-check my code and my reasoning ;-)
Execution is pretty fast (a blink of an eye on Free42, will run it on my 42S when I'm reasonably certain it is correct)

Cheers, Werner

[J-F Garnier]

(11-29-2022 08:52 AM)Werner Wrote:  I still have to double-check my code and my reasoning ;-)

At least you have a reasoning :-)
Clearly brute force is not the way to go, but I can't figure out a short-cut.
But knowing that you got a solution so quickly is already a clue.

J-F

[John Keith]

I wouldn't call f(n) an "unassuming series", it's really quite a complex and interesting series. Being recursive it can take some time to compute for larger numbers but there are some shortcuts we can use to speed up the computation.

First of all, d(n), the number of bits needed to represent n, is a very simple sequence that can be computed rapidly with a simple loop. Secondly, there is an underlying fractal pattern in f(n) and to compute k terms one need only compute approximately 2 * LOG2(k) terms recursively, and the rest can be filled in by a simple loop needing only 1 addition per term.

This is still basically brute force and I don't have a complete program ready yet but at least the germ of an idea.

[pinkman]

Are digressions allowed in the main thread?
Far from an answer, but I had a good surprise seeing the Prime (which will be consider as a 40 yo advanced calculator in 2053) can handle recursive functions in the apparently simple but yet powerful “Define” feature.
   

[ThomasF]

(11-29-2022 08:52 AM)Werner Wrote:  Well, I have a result..

2.08637766501

Hmm ... I also thought I had a result. I did some brute force and found a way to simplify, but maybe I'm headed somewhere else on a bad track ...
The track I followed led me closer to 1.96643... and should be possible to fit in a RPN version.

But, I just have to check some more to see if I'm on the wrong path this time

Cheers,
Thomas

[Werner]

Ran my program on my 42S, and produced the exact same result in 55 seconds.
Werner

[J-F Garnier]

Since we are several to have tried the brute force approach, here are my results that show that it's a no go.

I rewrote 1/f(n) to 1/(n.f(d)) and used the fact that all the n values between 2^d and 2^(d+1)-1 have the same number of binary digits d+1 so the sum of these values can be factored with the same f(d+1).
Translated into code (sorry for the change of notation):
J1=2^(K-1) @ J2=2*J1-1
S=0 @ FOR J=J1 TO J2 @ S=S+1/J @ NEXT J
S1=S1+S/F(K)

I also noticed that the intermediate value S converges to LOG(2), so to have an idea of how the series behaves I replace the partial sum S with LOG(2) for K>=15, i.e. for the f(n) with n>2^15.
I went up to K=4000, i.e. numbers with 4000 binary digits (> 10^1200) , and the sum has still not converged.
I can just say that the sum is larger than 1.95, which is compatible with both Werner and Thomas prelim results.

Any confirmation of my analysis is welcome!

10 ! SRC12C
20 L=4000
30 DIM F(L)
40 ! calculate the first L F(I)
50 F(1)=1 @ F(2)=2
60 FOR I=3 TO L
70 D=INT(LOG(I)/LOG(2))+1
90 F(I)=I*F(D)
100 NEXT I
120 !
130 S1=1/F(1)+1/F(2)+1/F(3)
140 FOR K=3 TO L
142 S=LOG(2) ! approx value used for K>=15
145 IF K>=15 THEN 210
150 J1=2^(K-1) @ J2=2*J1-1
170 S=0 @ FOR J=J1 TO J2 @ S=S+1/J @ NEXT J
210 S1=S1+S/F(K)
220 IF K<15 OR MOD(K,500)=0 THEN DISP K;S;S1
230 NEXT K
235 END

>RUN
K, partial sum S, sum S1
3 .759523809524 1.79325396826
4 .725371850372 1.82347779536
5 .709016202207 1.84711166877
6 .701020708269 1.86658446622
7 .697068688885 1.88318133976
8 .695104120223 1.88680167372
9 .694124696722 1.89001521398
10 .693635700225 1.89290536273
11 .693391380792 1.89553184523
12 .693269265777 1.89793903018
13 .69320821942 1.9001608514
14 .693177699089 1.90222388027
500 .69314718056 1.95019974933 (S estimated from here)
1000 .69314718056 1.95220716995
1500 .69314718056 1.95327727197
2000 .69314718056 1.95403237901
2500 .69314718056 1.95457439294
3000 .69314718056 1.9550131171
3500 .69314718056 1.95538406366
4000 .69314718056 1.95570539887

J-F

[ThomasF]

(11-30-2022 09:15 AM)J-F Garnier Wrote:  Any confirmation of my analysis is welcome!

Hi J-F,

I can't confirm, but that was the same track I was using.

Realizing that f(d) is constant between 2^d and 2^(d+1)-1, and if moved out of the summation we are left with a harmonic series, i.e we could just sum over something like this: 1/f(d) * (H(2^(d+1)-1) - H(2^d-1)) with d=1..K, and testing that hypothesis converges to something close to 1.96643...

At least we get the similar result - is it the right track ... ? Maybe

Edit: Corrected the formula - got one reciprocal too many ...

Cheers,
Thomas

[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

[PeterP]

Nicely done, Werner!

I got to the point where I had the sum on the left and the ln(2) times the partial sum on the right but I could not figure out the trick of adding the sum starting from 1 - K-1 back on the right side. Neat!

I used the in(ln2(x))+1 formula for the number of digits and a different approximation for Hn and but it has similar limits for accuracy. With a few tweaks your amazing code would fit the 41 as well, which I was shooting for.

[Fernando del Rey]

Wow! I've been working on-and-off on this new problem, and can't get beyond the useless brute force approach.

I'm not being able to find a good trick to speed up convergence. I've tried to transform the series into an equivalent series with faster convergence, but so far I did not get anywhere.

I have not given up, however, I fear that this time I won't be among the winners circle

[Valentin Albillo]

.
Hi, FLISZT, Werner, J-F Garnier, John Keith and Fernando del Rey,

A few assorted comments to what all of you said in your posts:

FLISZT Wrote:Yes I'm ready. With the 9th Symphony of Bruckner broadcasted at the radio / Eugen Jochum / Bavarian Radio Orchestra (1954) and now the Concerto for Piano & Orch #3 / Robert Casadesus, New York Philh, Dimitri Mitropoulos (live, 1957)… how could it not be the case?!

Hehe, another classical music lover, like me. There was a time many years ago when I listened to classical music exclusively, nothing else appealed to me; it lasted for a few years, then I diversified again. By the way, I guess that your FLISZT alias honors Franz Liszt, amirite ?

FLISZT Wrote:Here is my attempt:

Very nice, thank you for going ahead and posting it here. But that is just the definition of the recursive f(n), you still need to provide the code which computes the sum ... ... and btw, thanks for editing that code panel out.

Werner Wrote:Well, I have a result.. [...] I still have to double-check my code and my reasoning ;-) Execution is pretty fast (a blink of an eye on Free42, will run it on my 42S when I'm reasonably certain it is correct)

Ok, please post your code when you see fit but do try to avoid spoiling the fun for other people working on their own solutions right now. Oh, I see you've posted it already ...

Werner Wrote:Time for my reasoning and code ;-) [...] To calculate the H(n) [...] (thanks Albert!)

Thanks a lot for your very detailed reasoning and nice HP-42S RPN code, Werner, your teaming with Albert (Chan, I suppose, who avoids posting here for whatever reason) has indeed delivered the goods: code, result, timing.

The only thing remaining is to ascertain that your result is correct, which will be most likely the case if someone else reproduces it using a wholly different program. We'll see ...

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.

John Keith Wrote:I wouldn't call f(n) an "unassuming series", [...]

(tongue in cheek) "Write a sumptuous program to summarily sum this unassuming yet serious series:"

John Keith Wrote:it's really quite a complex and interesting series. Being recursive it can take some time to compute for larger numbers [...]

Not really, unless the numbers are humongously large. For instance, f(1,000,000,000) needs only the initial call and 4 recursive calls, if I counted'em right.

John Keith Wrote:This is still basically brute force and I don't have a complete program ready yet but at least the germ of an idea.

Good. Do persevere and post your completed program when ready, please. Remember: post code (no panels), sample run, result, timing and comments if at all possible.

Fernando del Rey Wrote:I've been working on-and-off on this new problem, and can't get beyond the useless brute force approach. [...] I've tried to transform the series into an equivalent series with faster convergence, but so far I did not get anywhere. I have not given up, however[...]

Don't despair, Fernando, it's just a matter of perseverance and lo and behold, the correct "trick" comes to mind in a flash. There's still plenty of time for you to produce a correct solution before I post my original one.


Well, thanks to all of you for your interest, you still have a number of days before I post my 6-line original solution.

Best regards.
V.

[John Keith]

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

John Keith Wrote:it's really quite a complex and interesting series. Being recursive it can take some time to compute for larger numbers [...]

Not really, unless the numbers are humongously large. For instance, f(1,000,000,000) needs only the initial call and 4 recursive calls, if I counted'em right.

John Keith Wrote:This is still basically brute force and I don't have a complete program ready yet but at least the germ of an idea.

Good. Do persevere and post your completed program when ready, please. Remember: post code (no panels), sample run, result, timing and comments if at all possible.


Best regards.
V.

Unfortunately my program was just a brute force one and as many have found, brute force is a dead end for this challenge. After summing 2^17 terms and seeing no convergence, it was obvious that I was getting nowhere. 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...

[FLISZT]

Hi All!

(11-30-2022 07:01 PM)Valentin Albillo Wrote:  Hehe, another classical music lover, like me. There was a time many years ago when I listened to classical music exclusively, nothing else appealed to me; it lasted for a few years, then I diversified again. By the way, I guess that your FLISZT alias honors Franz Liszt, amirite ?

Yes, rather a "mélomane" in general and a classical music lover in particular. I have two old souvenirs with classical music but out of topic. I have listened to different kinds of music. But unlike you, for the last few years I have been listening almost exclusively to classical music (symphonies, symphonic poems... and pieces created for the piano.)
Indeed, I like Liszt, but I could have choosen an other name like "Saint-Saëns" which sounds very "français" (so many composers created masterpieces) or "Franz Lisp" (I read that this language did exist! ) to be more connected to the forum.

(11-30-2022 07:01 PM)Valentin Albillo Wrote:  But that is just the definition of the recursive f(n)

Yes, I realized my mistake later… and so I understood why I was the 1st to publish a code when there are so many math and programming gurus in this forum!
A bit like thinking to be ahead in a F1 Grand-Prix with an old 2CV… but already being almost one lap late, 1 min after departure.
When I "draw" a calculator, it's basically only to check my accounts or to try this kind of "game". So I'm going to search with the "help" of my old and few skills in math.

(11-30-2022 07:01 PM)Valentin Albillo Wrote:  .. and btw, thanks for editing that code panel out.

Don't mention it!
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.
On the other hand, the layout is much more adapted to a structured language like RPL. Allowing small code panels (with no scrolling) could be a solution… or not. Just a thought.

[Werner]

(11-30-2022 07:01 PM)Valentin Albillo Wrote:  Ok, please post your code when you see fit but do try to avoid spoiling the fun for other people working on their own solutions right now. Oh, I see you've posted it already ...

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

Quote:Thanks a lot for your very detailed reasoning and nice HP-42S RPN code, Werner, your teaming with Albert has indeed delivered the goods: code, result, timing.

I didn't team up with Albert, I just used his formula for H(n), which is shorter and more accurate than the one you find elsewhere, and which he seems to be able to produce off the cuff. (Though in this case, as he pointed out to me in the meantime, using the original formula with a few more terms would've been better as it contains LN(2), which would then be cancelled out)

Quote:The only thing remaining is to ascertain that your result is correct, which will be most likely the case if someone else reproduces it using a wholly different program. We'll see ...

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

Cheers, Werner

[Werner]

Using the original asymptotic series H(n) = ln(n) + gamma + 1/(2*n) - 1/(12*n^2) + 1/(120*n^4) - ... results in:
H(n-1)-H(n/2-1) - ln(2) = 1/(2n) + 1/(4n^2) - 1/(8n^4) + 1/(4n^6) - 17/(16n^8) + .. (thanks^2, Albert!).
This formula is so accurate that my stopping criterion needed to be changed to 1+x = 1 ;-) But 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.

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

Cheers, Werner