Programming exercise (RPL/RPN) - Reciprocal Fibonacci Constant

[Gerson W. Barbosa]

Quoting from Wikipedia:

"The reciprocal Fibonacci constant, or ψ, is defined as the sum of the reciprocals of the Fibonacci numbers:

\(\psi = \sum_{k=1}^{\infty} \frac{1}{F_k} = \frac{1}{1} + \frac{1}{1} + \frac{1}{2} + \frac{1}{3} + \frac{1}{5} + \frac{1}{8} + \frac{1}{13} + \frac{1}{21} + \cdots.\) "

Our task is to write a program, the shortest the best, to compute the partial sums of this series from k=1 up to a given n. For instance, on the HP 50g, assuming the program is named RFC:

1 RFC --> 1.
2 RFC --> 2.
3 RFC --> 2.5
4 RFC --> 2.83333333333
5 RFC --> 3.03333333333
6 RFC --> 3.15833333333
7 RFC --> 3.23525641025

Convergence to d-digit results occurs when n is around ⌈(d*ln(100) - ln(20))/(2*ln(φ)⌉, where φ is the golden ratio (1.61803398875...). Thus, on the HP-41, we will need at least 46 terms for the exact 10-figure result:

45 XEQ ALPHA RFC ALPHA --> 3.359885665
46 XEQ ALPHA RFC ALPHA --> 3.359885666

On Free42, we can get at least 33 correct digits:

160 XEQ RFC --> 3.35988566624317755317201130291892(3)

By the way, this might be a breeze on the wp34s, which has FIB built in :-)

As a reference, my counts are

HP 50g: 50 bytes

HP-48G: 52.5 bytes

HP-42S: 28 bytes (HP-41 compatible)

HP-41CV: 33 bytes

These are only second (HP-41 & 42), or third (HP-48 and 50 g) attempts, so there surely is room for improvement.

Have fun!

Gerson.

[Paul Dale]

(02-16-2017 08:29 PM)Gerson W. Barbosa Wrote:  By the way, this might be a breeze on the wp34s, which has FIB built in :-)

And a summation command

[pre]LBL A
[Sigma] 00
RTN
LBL 00
FIB
1/x
[/pre]

I've not tried this code.


Pauli

[Gerson W. Barbosa]

(02-16-2017 10:24 PM)Paul Dale Wrote:  

(02-16-2017 08:29 PM)Gerson W. Barbosa Wrote:  By the way, this might be a breeze on the wp34s, which has FIB built in :-)

And a summation command

[pre]LBL A
[Sigma] 00
RTN
LBL 00
FIB
1/x
[/pre]

I've not tried this code.


Pauli

Missing only three instructions between LBL A and Sigma 00:

#001
SDR 003
+

I had to check Walter's Blue Book ( page 124 ).

160 A gives 34 correct digits in DOUBLE mode, BTW.

Gerson.

PS - Nine steps only! I had to use fourteen on the HP-41/42S.

[Ángel Martin]

(02-16-2017 11:44 PM)Gerson W. Barbosa Wrote:  PS - Nine steps only! I had to use fourteen on the HP-41/42S.

Will be down to only one step with the MCODE function "s1/FIB" ;-)
In preparation, thanks for the week-end exercise...

[Paul Dale]

(02-16-2017 11:44 PM)Gerson W. Barbosa Wrote:  PS - Nine steps only! I had to use fourteen on the HP-41/42S.

Thanks for the correction

I'm sure the fourteen step version will be a lot faster on the 34S than mine: FIB is a fairly expensive function and Sigma is a key stroke program that Kahan sums the terms.


Pauli

[Paul Dale]

How's this:

Code:

    STO 00
    CLR STK
    ISG X
    LBL 00
        STO+ Y
        1/x
        STO+ Z
        1/x
        x<>y
        DSE 00
            GTO 00
    RCL Z

Twelve steps, eighteen bytes. This doesn't use any functions they don't have, although the accuracy won't be so great.

In thirteen steps, twenty bytes and no register:

Code:


    0
    ENTER
    ENTER
    ISG X
    LBL 00
        STO+ Y
        1/x
        STO+ Z
        1/x
        x<>y
        DSE T
            GTO 00
    RCL Z

Pauli

[Werner]

How can the 42S code be 28 bytes, 41-compatible, yet the 41CV code is 33 bytes?

I have 27 bytes (excluding END, but including the three-letter label) on the 42 and 25 on the 41, stack-only. Of course.

Code:

>LBL "RFB"
 0
 STO ST Z
 1
>LBL 02
 1/X
 STO+ ST T
 X<> ST L
 STO+ ST Y
 X<>Y
 DSE ST Z
 GTO 02
 R^
 END

Cheers, Werner

[Guenter Schink]

(02-17-2017 08:03 AM)Werner Wrote:  How can the 42S code be 28 bytes, 41-compatible, yet the 41CV code is 33 bytes?

I have 27 bytes (excluding END, but including the three-letter label) on the 42 and 25 on the 41, stack-only. Of course.

Code:

>LBL "RFB"
 0
 STO ST Z
 1
>LBL 02
 1/X
 STO+ ST T
 X<> ST L
 STO+ ST Y
 X<>Y
 DSE ST Z
 GTO 02
 R^
 END

Günter

Cheers, Werner

Why excluding "END"? My program for the 42S is pretty close to yours. Also 27 bytes, including "END" and stack only, of course

Code:

01 LBL "RFC"
02 0
03 0
04 1
05 LBL 00
06 1/x
07 STO+ ST Z
08 1/X
09 STO+ ST Y
10 x<>Y
11 DSE ST T
12 GTO 00
13 RCL ST Z
14 END

Günter

[Werner]

I don't count the END because the 42S doesn't, either.

Code:

06 1/X
07 STO+ ST Z
08 1/X

Well, I would never do this ;-) Even if it doesn't seem to make any difference here.
Werner

[Guenter Schink]

(02-17-2017 12:15 PM)Werner Wrote:  I don't count the END because the 42S doesn't, either.

Code:

06 1/X
07 STO+ ST Z
08 1/X

Well, I would never do this ;-) Even if it doesn't seem to make any difference here.
Werner

Ups, you're right. That could give unexpected results. Replace by LASTx.

Günter

[Gerson W. Barbosa]

Here are mine, all stack-only of course, mainly because these run faster. However, numbered registers and local variable versions might be better if lower byte count is the only goal. Perhaps we can divide the solutions into two groups. So far for the HP-42S we have a few stack-only 27-byte solutions and one 25-byte solution, the latter using one number register.

Congratulations to all!

Gerson.


HP-50g (50 bytes)

« 0. 1. DUP2 5. ROLL

START SWAP OVER + ROT OVER INV + UNROT

NEXT DROP2

»



HP-48G (52.5 bytes)

« 0 1 DUP2 ROT 5 ROLL

START OVER + ROT
OVER INV + SWAP ROT

NEXT DROP2

»


HP-42S

00 { 28-Byte Prgm }
01>LBL "RFC"
02 0
03 0
04 1
05>LBL 00
06 +
07 LASTX
08 1/X
09 STO+ ST Z
10 X<> ST L
11 X<>Y
12 DSE ST T
13 GTO 00
14 RCL ST Z
15 .END.


HP-41

01>LBL 'RFC
02 0
03 0
04 1
05>LBL 00
06 +
07 LASTX
08 1/X
09 STO+ Z
10 X<> L
11 X<>Y
12 DSE T
13 GTO 00
14 RCL Z
15 .END.

[Gerson W. Barbosa]

(02-17-2017 08:03 AM)Werner Wrote:  How can the 42S code be 28 bytes, 41-compatible, yet the 41CV code is 33 bytes?

I think it's strange too, but that's what I got on my HP-41CV with CAT 1 and the printer on, after PACKING. Perhaps I should do a manual byte count.

Cheers,

Gerson.

PS

(02-17-2017 08:03 AM)Werner Wrote:  

Code:

>LBL "RFB"
 0
 STO ST Z
 1
>LBL 02
 1/X
 STO+ ST T
 X<> ST L
 STO+ ST Y
 X<>Y
 DSE ST Z
 GTO 02
 R^
 END

55 XEQ RFB --> 3.3598856622 (7.9 seconds)

[Gerson W. Barbosa]

(02-17-2017 07:09 AM)Paul Dale Wrote:  I'm sure the fourteen step version will be a lot faster on the 34S than mine: FIB is a fairly expensive function and Sigma is a key stroke program that Kahan sums the terms.

Indeed, 21+ seconds versus slightly less than one second. Only 33 correct digits, however.

Gerson.

[Gerson W. Barbosa]

(02-17-2017 07:29 AM)Paul Dale Wrote:  How's this:

Code:

    STO 00
    CLR STK
    ISG X
    LBL 00
        STO+ Y
        1/x
        STO+ Z
        1/x
        x<>y
        DSE 00
            GTO 00
    RCL Z

Twelve steps, eighteen bytes. This doesn't use any functions they don't have, although the accuracy won't be so great.

In thirteen steps, twenty bytes and no register:

Code:


    0
    ENTER
    ENTER
    ISG X
    LBL 00
        STO+ Y
        1/x
        STO+ Z
        1/x
        x<>y
        DSE T
            GTO 00
    RCL Z

Pauli

One more step, considering the initial LBL. 25 and 27 bytes on the HP-42S, respectively. 55 iterations in 9.0 s and 8.1 s, respectively.

Gerson.

[Claudio L.]

(02-17-2017 01:58 PM)Gerson W. Barbosa Wrote:  [font=Courier]HP-50g (50 bytes)

« 0. 1. DUP2 5. ROLL

START SWAP OVER + ROT OVER INV + UNROT

NEXT DROP2

»

Here's an alternative, using local variables and lists for readability:

Code:


« → N 
    «
     1 1
     3 N START DUP2 + NEXT
     N →LIST 
     INV ∑LIST
   »
»

I didn't measure on stock firmware since I have my calc with newRPL, but for comparison your code was 72 bytes in newRPL and the alternative is 80 bytes.
One glitch: it only works for N=3 and higher.
Just a different point of view.

[Thomas Okken]

(02-17-2017 02:09 PM)Gerson W. Barbosa Wrote:  

(02-17-2017 08:03 AM)Werner Wrote:  How can the 42S code be 28 bytes, 41-compatible, yet the 41CV code is 33 bytes?

I think it's strange too, but that's what I got on my HP-41CV with CAT 1 and the printer on, after PACKING. Perhaps I should do a manual byte count.

The HP-42S' byte count doesn't include the END, which of course is cheating, because that END (or .END.) really does take up three bytes.
That still doesn't explain the discrepancy, though, and there is another factor that should make the 41 version one byte shorter than the 42S version, not two bytes longer: on the 42S, numbers are always preceded by a null byte, while the 41 only inserts a null byte between consecutive number lines.

[Gerson W. Barbosa]

(02-17-2017 09:50 PM)Claudio L. Wrote:  

(02-17-2017 01:58 PM)Gerson W. Barbosa Wrote:  HP-50g (50 bytes)

« 0. 1. DUP2 5. ROLL

START SWAP OVER + ROT OVER INV + UNROT

NEXT DROP2

»

Here's an alternative, using local variables and lists for readability:

Code:


« → N 
    «
     1 1
     3 N START DUP2 + NEXT
     N →LIST 
     INV ∑LIST
   »
»

I didn't measure on stock firmware since I have my calc with newRPL, but for comparison your code was 72 bytes in newRPL and the alternative is 80 bytes.
One glitch: it only works for N=3 and higher.

59 bytes in RPL.

(02-17-2017 09:50 PM)Claudio L. Wrote:  Just a different point of view.

Also, a great idea. Thank you very much for your contribution!

This modified stack-only version is 1.5 bytes longer. The glitch still persists for N=1, but only because ∑LIST doesn't accept one-element lists, at least here in old RPN.

Gerson.

Code:


« 0. { } 1. 1. 5. ROLL
  START PICK3 + DUP UNROT + ROT
  NEXT ROT DROP2 INV ∑LIST
»

[Joe Horn]

(02-16-2017 08:29 PM)Gerson W. Barbosa Wrote:  7 RFC --> 3.23525641025

Many of the programs submitted in this thread suffer from avoidable roundoff error due to summing the reciprocals starting with the largest and ending with the smallest, which is always a Bad Thing to do. For example, the above value is incorrect (the last digit should be 6), but if the reciprocals are summed in the opposite order, beginning with 1/F(7) and ending with 1/F(1), then the correct result is obtained.

In general, the following RFC(n) when obtained by summing 1/F(1) through 1/F(n) are incorrect (in a 12-digit mantissa machine), but when obtained by summing 1/F(n) through 1/F(1) they are correct (rounded to 12 significant digits of course):

RFC 7; 8; 13-24; 27-29; 31-36; and all n>=38.

Both summation orders fail only for RFC 25 & 37; higher precision methods are required for those two values of n. For reference, here are their correct rounded values:
RFC(25) = 3.35986409965
RFC(37) = 3.35988559927

[Gerson W. Barbosa]

(02-17-2017 10:36 PM)Thomas Okken Wrote:  

(02-17-2017 02:09 PM)Gerson W. Barbosa Wrote:  I think it's strange too, but that's what I got on my HP-41CV with CAT 1 and the printer on, after PACKING. Perhaps I should do a manual byte count.

The HP-42S' byte count doesn't include the END, which of course is cheating, because that END (or .END.) really does take up three bytes.
That still doesn't explain the discrepancy, though, and there is another factor that should make the 41 version one byte shorter than the 42S version, not two bytes longer: on the 42S, numbers are always preceded by a null byte, while the 41 only inserts a null byte between consecutive number lines.

Thanks for your detailed explanation. I think I've found out what happened. I keyed the same program into my HP-41C and got 32 bytes, which is still the wrong byte count. In both calculators that was the last program in the catalog list. After entering another small program into the HP-41C I finally obtained 30 bytes. Apparently the CAT 1 method doesn't work for the most recent program in the calculator. I've been using an infrared module and an HP-82240B printer.

Gerson.

[Gerson W. Barbosa]

(02-18-2017 12:53 AM)Joe Horn Wrote:  

(02-16-2017 08:29 PM)Gerson W. Barbosa Wrote:  7 RFC --> 3.23525641025

Many of the programs submitted in this thread suffer from avoidable roundoff error due to summing the reciprocals starting with the largest and ending with the smallest, which is always a Bad Thing to do. For example, the above value is incorrect (the last digit should be 6), but if the reciprocals are summed in the opposite order, beginning with 1/F(7) and ending with 1/F(1), then the correct result is obtained.

In general, the following RFC(n) when obtained by summing 1/F(1) through 1/F(n) are incorrect (in a 12-digit mantissa machine), but when obtained by summing 1/F(n) through 1/F(1) they are correct (rounded to 12 significant digits of course):

RFC 7; 8; 13-24; 27-29; 31-36; and all n>=38.

Both summation orders fail only for RFC 25 & 37; higher precision methods are required for those two values of n. For reference, here are their correct rounded values:
RFC(25) = 3.35986409965
RFC(37) = 3.35988559927

Yes, I do remember this fact from another program exercise I first did in 1984 when I was studying Pascal, exercise 8.3 in Systematic Programming, An Introduction, by Prof. Niklaus Wirth. :-)

That's one of the reasons, along with Kahan's sum in the built-in summation command, that Paul Dale's first wp34s program gives all 34 correct digits. I had tried REVLIST after INV in Claudio L.'s program and noticed some improvement in accuracy.

Doing the sum backwards make programs longer, however. But your are right, accuracy-wise, that's the method that should be used.

Thank you!

Gerson.