Mardi Gras True Fibs

[Gerson W. Barbosa]

This is a follow-up of the recent thread about the Reciprocal Fibonacci Constant.

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.\) "

That was an exercise on computing this series up to a certain number of terms. A number of you have contributed interesting codes and ideas, to which I am grateful.

It has been mentioned that the computation of the constant to d digits requires at least ⌈(d*ln(100) - ln(20))/(2*ln(φ)⌉ terms, where φ is the golden ratio ( (1 + √5)/2 ), using the plain series. Thus, for 100 digits we would need to sum up to 476 terms of the series. The demonstration is trivial and is left as an exercise for the reader.
476 terms to compute 100 digits looks somewhat inefficient, so my goal has shifted to finding a way to lower this number. Here I share what I have found. Again, without a proof, but only because I have none:

\[\psi = \frac{1}{F_{1}}+\frac{1}{F_{2}}+\frac{1}{F_{3}}+\cdots +\frac{1}{F_{n-1}}+\frac{1}{F_{n}}+\frac{1}{F_{n-1}-\frac{F_{1}^{2}}{F_{n+2}\cdot F_{1}+F_{n-1}\cdot F_{2}-\frac{F_{2}^{2}}{F_{n+2}\cdot F_{2}+F_{n-1}\cdot F_{3}-\frac{F_{3}^{2}}{F_{n+2}\cdot F_{3}+F_{n-1}\cdot F_{4}-\frac{F_{4}^{2}}{F_{n+2}\cdot F_{4}+F_{n-1}\cdot F_{5}-\frac{F_{5}^{2}}{F_{n+2}\cdot F_{5}+F_{n-1}\cdot F_{6}-\frac{F_{6}^{2}}{F_{n+2}\cdot F_{6}+F_{n-1}\cdot F_{7}-... }}}}}}}\]

Code:


%%HP: T(3)A(R)F(.);
\<< 0 R\<-\->F 1 R\<-\->F 1 4 PICK
  START DUP 4 ROLLD DUP ROT FADD
  NEXT FADD 4 PICK ROT 2 + ROLLD LASTARG ROLLD LASTARG 2 - 0 1 ROT
  START SWAP FINV FADD
  NEXT 4 ROLLD 1 4 PICK
  START DUP 4 ROLLD DUP ROT FADD
  NEXT DROP2 1 - 1
  FOR i i FIBn R\<-\->F FSQ SWAP FDIV FNEG FADD -1
  STEP FINV FADD ZZ\<-\->F DROP \->STR DUP HEAD "." + SWAP TAIL + DUP SIZE 1 - "  " REPL
\>>

FIBn:

%%HP: T(3)A(R)F(.);
\<< 5. \v/ DUP 1. + 2. / ROT ^ SWAP / .5 + FLOOR
\>>

This HP 50g user RPL program, using the LongFloat library, requires two parameters. The first is the number of terms of the continued fraction and the second is n in the expression above, the number of terms of the original series.

The following is an evaluation of the expression above with various parameters. When we compare the results with the true reciprocal Fibonacci constant in the first line, we notice that best results are achieved with approximately equal numbers of terms of the continued fraction and terms of the plain series. In this case we have obtained 99 correct decimals using only 30 terms in the last line. 30 is about 6.3% of 476. Not bad!

       3.359885666243177553172011302918927179688905133731968486495553815325130318996683​383615416216456790087

  2 12 3.3598856662421299884649892392688576236687368131807676633574000702972674942939991993079972711​21544920

  2 20 3.359885666243177553167433281787748928648776577928290964879706040408646726702491028523646222602832650​

  2 30 3.359885666243177553172011302918764514065619978389958675293504780890783041151282869607821301725790765

  2 40 3.359885666243177553172011302918927179688899353919213085077297503281892745245057092283592271334021082

  2 50 3.359885666243177553172011302918927179688905133731968281128034773312024205107370483101951938492442443

  4 10 3.3598856662431775529330373071740194358904425219090280535498309533751379799543475599629279294234996​62

  4 20 3.359885666243177553172011302918927179651795581225695045861072258099578065882672826328709609425804709

  4 30 3.359885666243177553172011302918927179688905133731968486489791485719247471989816805452564095340121083

  4 40 3.359885666243177553172011302918927179688905133731968486495553815325130318995788615553566782514358382

  4 50 3.359885666243177553172011302918927179688905133731968486495553815325130318996683​383615416216456651148

 10 20 3.359885666243177553172011302918927179688905133731968486495553815325130318996683​383615416201795067149

 11 20 3.359885666243177553172011302918927179688905133731968486495553815325130318996683​383615416216456790047

 14 14 3.359885666243177553172011302918927179688905133731968486495553815325130318996683​383615414413089430501

 15 14 3.359885666243177553172011302918927179688905133731968486495553815325130318996683​383615416216455221736

 14 16 3.359885666243177553172011302918927179688905133731968486495553815325130318996683​383615416216456790085


Hopefully there are not HP-49/50g-exclusive commands in the following RPL program. It is more difficult to be used with the LongFloat library, though.

Code:


%%HP: T(3)A(R)F(.);
\<< 0. 1. 1. 4. PICK
  START DUP 4. ROLLD DUP ROT +
  NEXT + 4. PICK ROT 2. + ROLLD LASTARG ROLLD LASTARG 2. - \->LIST REVLIST INV \GSLIST 4. ROLLD 1. 4. PICK
  START DUP 4. ROLLD DUP ROT +
  NEXT DROP2 \->LIST REVLIST DUP SIZE 1. - ROT SWAP 0. 1. 1. 4. PICK
  START DUP 4. ROLLD DUP ROT +
  NEXT DROP2 \->LIST REVLIST SQ NEG 0. + ROT 0. + OVER SIZE 1. DUP ROT
  START GETI SWAP 1. - 4. ROLL SWAP GETI 4. ROLL / 4. ROLL ROT GETI SWAP 1. - SWAP 4. ROLL + PUTI 1. -
  NEXT DROP DUP SIZE 1. - GET INV SWAP DROP +
\>>

Comments, corrections and code optimizations are welcome.

Definitely it looks like I have to find something more interesting to do during Easter and Carnival holidays ( Mardi Gras Basic Trigs (HP-12C), Easter Sunday Basic Trigs (HP-12C), Easter Sunday Trigs ( rpn38-CX) ).

Next year I think I'll go to Rio :-)

Edited to fix a couple of typos and minor mistakes.

[Don Shepherd]

Gerson, you continue to amaze!

[Gerson W. Barbosa]

I forgot to mention that n (the number of terms of the original series) has to be even. The continued fraction has to have at least two terms. The latter is a limitation of the program only. Notice the continuous fraction terms involved one division and one subtraction only. Let's take a look at a numerical example:


\[\frac{1}{1}+\frac{1}{1}+\frac{1}{1-\frac{1^{2}}{4-\frac{1^{2}}{5-\frac{2^{2}}{9-\frac{3^{2}}{14-\frac{5^{2}}{23-\frac{8^{2}}{37}}}}}}}=3.35988549037\]

That's the result we get when running either program with arguments 7 and 2. In this case, n = 2 and F(n-1) = F1 = 2; F(n+2)*F(1) + F(n-1)*F(2) = F(4)*F(1) + F(1)*F(2) = 3*1 + 1*1 = 4. Once these have been calculated the next terms are obtained by simple addition: 5, 9, 14, 23, 37..., that is a new Fibonacci sequence with initial terms 1 and 4. Notice the first program uses FIBn only because this made the use of the LongFloat library easier.

Using an emulator 1000 digits can be computed in less than 10 seconds @ 2.6 GHz, running the program with parameters 49 and 48. This is only a small fraction of the terms required by the original series, 97/4800 (about 2%). Or about 20% of the somewhat more complicated terms used by the best solution here.

3.359885666243177553172011302918927179688905133731968486495553815325130318996683​38361541621645679008729704534292885391330413678901710088367959135173307711907858​03335503325077531875998504871797778970060395645092153758927752656733540240331694​41799293934610992626257964647651868659449710216558984360881472693249591079473873​67337852332687749976272775794685367691854198146766874299876738209691390121772202​44052081510942649349513745416672789553444707777758478025963407690748474155579104​20067501520341070533528512979263524206226753756805576195566972084884385440798332​42928513680708275226625797511886464640967374615723872362955620536122030246354092​52678424224347036310363201466298040249015578724456176000319551987905969942029178​86694917480809674652368265408693839906987321175216695706385941181455364736426878​24629261666501000989038048233595198931461501082887263928876699171493040530577455​74321561167298985617729731395370735291966884327898022165047585028091806291002444​27701746024104041778606919006503714283529

All digits correct according to that reference.

PS -

I've noticed results don't change when the parameters are swapped. Also, since apparently the maximum efficiency is reached when they're the same, there is no need to compute two different Fibonacci sequences with different sizes as these test programs have been doing. This will improve efficiency even more.

[Gerson W. Barbosa]

Here is a UserRPL solution that works on the HP 49G/50g and doesn't rely on external arbitrary precision libraries:

Code:


%%HP: T(3)A(R)F(.);
DIR
  RFC
  \<< PUSH RAD -105 CF -3 CF 0 1 1 4 PICK
    START DUP 4 ROLLD DUP ROT +
    NEXT + 4 PICK ROT 2 + ROLLD LASTARG ROLLD LASTARG 2 - \->LIST REVLIST DUP 4 ROLLD DUP SIZE 4 ROLLD INV \GSLIST 4 ROLLD 1 4 PICK
    START DUP 4 ROLLD DUP ROT +
    NEXT DROP2 \->LIST REVLIST SWAP ROT TAIL SQ NEG 0 + ROT 0 + OVER SIZE 1 DUP ROT
    START GETI SWAP 1 - 4 ROLL SWAP GETI 4 ROLL / 4 ROLL ROT GETI SWAP 1 - SWAP 4 ROLL + PUTI 1 -
    NEXT DROP DUP SIZE 1 - GET INV NIP + EXPAND \->STR "/" " " SREPL DROP "'" "" SREPL DROP OBJ\-> OVER \->STR SIZE .89 ^ 1.209 * 1 - IP 
    R\->I ALOG ROT OVER * UNROT * DUP2 DUP SIZE R\->I ALOG SWAP - * SWAP IQUOT + \->STR DUP HEAD -51 FC? { "." } { "," } IFTE + SWAP TAIL + POP
  \>>
  d
  \<< DUP \v/ 1.5 * CEIL DUP 2 MOD + RFC 1 ROT 1 + SUB
  \>>
END

The following table and information might give an idea of the method.


 n  | d   |  N   |  D  |  k 
----+-----+------+-----+----
  8 |  30 |  15  |  14 |  15
 10 |  45 |  24  |  24 |  20
 12 |  64 |  36  |  35 |  28
 14 |  87 |  50  |  49 |  38
 16 | 113 |  64  |  64 |  48
 18 | 142 |  83  |  82 |  59
 20 | 174 | 100  | 100 |  73
 22 | 211 | 123  | 123 |  87
 40 | 685 | 419  | 418 | 266
 48 | 981 | 612  | 611 | 369
 50 |1064 | 664  | 663 | 400


n: number of terms of the regular term = number of terms of the continued fraction
d: number of correct digits
N: number of digits of the numerator of the simplified fraction
D: number of digits of the denominator of the simplified fraction
k: exponent of 10 used in the conversion from simplified fraction to decimal fraction

N1 = N*10^k
D1 = D*10^k

N2 = (10_complement(D1)+1)*N1) idiv (D1)
D2 = 10^(2*k)

k ~ 1.2093*N^0.88996 ( in the program, 1.209*N^0.89 - 1 )
n ~ (2.7735*d^0.51461)/2 ( in the program, ceil[1.5*sqrt(d)] )

The first program evaluates the Reciprocal Fibonacci Constant for n terms of the regular series and n terms of the continued fraction. The second program evaluates it to d decimal digits, including the integer part. The latter occasionally will fail since a better curve-fitting is still required. For instance, 1001 d will return a 977-character string rather than the 1001 digits we would expect. Using ceil[1.5*sqrt(d)] + 2 would help in this particular case, but would be worse overall.

Examples:

 2  RFC  -->    3.3  (2  digits,  1.06  s)
 8  RFC  -->    3.35988566624317755317201130  (28  digits,  2.62  s)
16  RFC  -->    3.359885666243177553172011302918927179688905133731968486495553815325130318996683​383615416216456790087297045342928  (112  digits,  13.40  s)

  50  d  -->    3.3598856662431775531720113029189271796889051337319  (6.00  s)
 100  d  -->    3.359885666243177553172011302918927179688905133731968486495553815325130318996683​383615416216456790087  (13.53  s)


1500 digits under one hour, 1000 digits under 20 minutes. Times on the HP 50g (about three times longer on the HP 49G).

Edited. N and D are the number of digits of the numerator and the denominator of the resulting irreducible fraction, not the numerator and the denominator themselves, in case some might have wondered. Sorry for the lack of attention.

—————————-

Updated version using the FXND command:

« PUSH RAD -105 CF -3 CF DUP √ 1.5 * CEIL DUP 2 MOD + 0 1 1 4 PICK
START DUP 4 ROLLD DUP ROT +
NEXT + 4 PICK ROT 2 + ROLLD LASTARG ROLLD LASTARG 2 - →LIST REVLIST DUP 4 ROLLD DUP SIZE 4 ROLLD INV ∑LIST 4 ROLLD 1 4 PICK
START DUP 4 ROLLD DUP ROT +
NEXT DROP2 →LIST REVLIST SWAP ROT TAIL SQ NEG 0 + ROT 0 + OVER SIZE 1 DUP ROT
START GETI SWAP 1 - 4 ROLL SWAP GETI 4 ROLL / 4 ROLL ROT GETI SWAP 1 - SWAP 4 ROLL + PUTI 1 -
NEXT DROP DUP SIZE 1 - GET INV NIP + EXPAND FXND OVER →STR SIZE .89 ^ 1.209 * 1 - IP R→I ALOG ROT OVER * UNROT * DUP2 DUP SIZE R→I ALOG SWAP - * SWAP IQUOT + →STR DUP HEAD -51 FC? { "." } { "," } IFTE + SWAP TAIL + 1 ROT 1 + SUB POP
»

The argument is the number of figures, as in the previous examples.

[Gerson W. Barbosa]

This thread is more than one year old, but I have decided to updated it for the sake of completeness.

Firstly, the the following revised formula is not restricted to an even number of terms anymore:

\(\psi \approx \frac{1}{F_{1}}+\frac{1}{F_{2}}+\frac{1}{F_{3}}+\cdots +\frac{1}{F_{n-1}}+\frac{1}{F_{n}}+\frac{1}{F_{n-1}+\frac{s\cdot F_{1}^{2}}{F_{n+2}\cdot F_{1}+F_{n-1}\cdot F_{2}+\frac{s\cdot F_{2}^{2}}{F_{n+2}\cdot F_{2}+F_{n-1}\cdot F_{3}+\frac{s\cdot F_{3}^{2}}{F_{n+2}\cdot F_{3}+F_{n-1}\cdot F_{4}+\frac{s\cdot F_{4}^{2}}{F_{n+2}\cdot F_{4}+F_{n-1}\cdot F_{5}+\frac{\ddots }{\ddots F_{n+2}\cdot F_{n-1}+F_{n-1}\cdot F_{n}+\frac{s\cdot F_{n}^{2}}{F_{n+2}\cdot F_{n}+F_{n-1}\cdot F_{n+1} }}}}}}}\)

where \(s = \left \{ _{ 1, \mod(n,2)=1}^{-1, \mod(n,2)=0 } \right.\)

and

\(n = \left \lceil \sqrt{d}\left ( \varphi -\frac{1}{5\sqrt[8]{d}} \right )-\frac{1}{2} \right \rceil\)

where d is the number of correct decimal digits obtained from a given n.

and \(\varphi = \frac{1+\sqrt{5}}{2}\) is the golden ratio.

The formula \(n=f(d)\) above has been obtained empirically, but the fit is so close it might as well be exact:

      n           d       sqrt(d)*(phi - 1/(5*d^(1/8))) - 1/2   relative error (%)

     12         68.987                  11.96                         0.33
     16        120.326                  16.04                         0.27
     20        182.984                  19.98                         0.12
     22        220.489                  22.01                         0.06
     26        303.895                  26.00                         0.00
     30        400.909                  30.00                         0.01
     34        511.683                  34.03                         0.08
     38        635.227                  38.03                         0.08
     40        702.522                  40.05                         0.12
     44        846.221                  44.06                         0.14
     46        922.801                  46.06                         0.14
     48       1000.700                  48.02                         0.04

Secondly, I have provided a Decimal BASIC program. The algorithm is better than the one used in the previous RPL version, also the code might be easier to convert to other programming languages.

OPTION ARITHMETIC DECIMAL_HIGH
! Reciprocal Fibonacci Constant
DEF SQRT(x) =  EXP(LOG(x)/2)                                   ! standard precision SQR as
DEF Fib(x) = INT((phi^x - (-phi)^(-x))/SQRT(5) + 0.01)         ! full precision is not re-       
INPUT  PROMPT "Number of digits: ":f                           ! quired for low values of k                      
LET t = TIME                          
LET phi = (1 + SQRT(5))/2
LET k = CEIL(SQRT(f)*(phi - 1/(5*SQRT(SQRT(SQRT(f))))) - 1/2)
LET s = 2*MOD(k,2) - 1                                         ! if k is even then s = -1 else s = 1
LET sr = 0
LET a = Fib(k) 
LET b = Fib(k + 1)
LET cf = 0
LET u = Fib(k - 1)
LET v = Fib(1)*Fib(k + 2) + Fib(2)*Fib(k - 1)                  ! this expression evaluates correct-  
LET d = v*Fib(k) + u*Fib(k - 1)                                ! ly  for k < 70 using standard pre-   
LET e = v*Fib(k + 1) + u*Fib(k)                                ! cision  SQRT,  which would suffice  
FOR i = 1 TO k                                                 ! for 2000 digits.  Anyway,  Decimal  
   LET sr = sr + 1/a                                           ! BASIC is limited to 1000 digits 
   LET w = a
   LET a = b - a
   LET b = w
   LET n = s*b*b 
   LET cf = n/(cf + d)
   LET w = d
   LET d = e - d  
   LET e = w
NEXT i
LET cf = 1/(cf + d)
LET r = TIME - t                            
LET r$ = STR$(sr + cf)
PRINT
PRINT r$(1:2);
FOR i = 3 TO f + 1
   PRINT r$(i:i);
   IF MOD((i - 2),10) = 0 THEN PRINT " ";
   IF MOD((i - 2),50) = 0 THEN 
      PRINT
      PRINT "  ";
   END IF
NEXT i
IF MOD (i - 3,50) <> 0  OR f = 0 THEN PRINT
PRINT
PRINT "Runtime: ";
PRINT  USING "0.##": r;
PRINT " seconds"
END


Number of decimal places: 1001

3.3598856662 4317755317 2011302918 9271796889 0513373196 
  8486495553 8153251303 1899668338 3615416216 4567900872 
  9704534292 8853913304 1367890171 0088367959 1351733077 
  1190785803 3355033250 7753187599 8504871797 7789700603 
  9564509215 3758927752 6567335402 4033169441 7992939346 
  1099262625 7964647651 8686594497 1021655898 4360881472 
  6932495910 7947387367 3378523326 8774997627 2775794685 
  3676918541 9814676687 4299876738 2096913901 2177220244 
  0520815109 4264934951 3745416672 7895534447 0777775847 
  8025963407 6907484741 5557910420 0675015203 4107053352 
  8512979263 5242062267 5375680557 6195566972 0848843854 
  4079833242 9285136807 0827522662 5797511886 4646409673 
  7461572387 2362955620 5361220302 4635409252 6784242243 
  4703631036 3201466298 0402490155 7872445617 6000319551 
  9879059699 4202917886 6949174808 0967465236 8265408693 
  8399069873 2117521669 5706385941 1814553647 3642687824 
  6292616665 0100098903 8048233595 1989314615 0108288726 
  3928876699 1714930405 3057745574 3215611672 9898561772 
  9731395370 7352919668 8432789802 2165047585 0280918062 
  9100244427 7017460241 0404177860 6919006503 7142835295 
  
Runtime: 0.06 seconds

[Gerson W. Barbosa]

(06-14-2018 12:41 AM)Gerson W. Barbosa Wrote:  \(\psi \approx \frac{1}{F_{1}}+\frac{1}{F_{2}}+\frac{1}{F_{3}}+\cdots +\frac{1}{F_{n-1}}+\frac{1}{F_{n}}+\frac{1}{F_{n-1}+\frac{s\cdot F_{1}^{2}}{F_{n+2}\cdot F_{1}+F_{n-1}\cdot F_{2}+\frac{s\cdot F_{2}^{2}}{F_{n+2}\cdot F_{2}+F_{n-1}\cdot F_{3}+\frac{s\cdot F_{3}^{2}}{F_{n+2}\cdot F_{3}+F_{n-1}\cdot F_{4}+\frac{s\cdot F_{4}^{2}}{F_{n+2}\cdot F_{4}+F_{n-1}\cdot F_{5}+\frac{\ddots }{\ddots F_{n+2}\cdot F_{n-1}+F_{n-1}\cdot F_{n}+\frac{s\cdot F_{n}^{2}}{F_{n+2}\cdot F_{n}+F_{n-1}\cdot F_{n+1} }}}}}}}\)

where \(s = \left \{ _{ 1, \mod(n,2)=1}^{-1, \mod(n,2)=0 } \right.\)

Or, in a more compact way:

\(\psi \approx \frac{1}{F_{1}}+\frac{1}{F_{2}}+\frac{1}{F_{3}}+\cdots +\frac{1}{F_{n-1}}+\frac{1}{F_{n}}+\frac{1}{G_1+\frac{s\cdot F_{1}^{2}}{G_2+\frac{s\cdot F_{2}^{2}}{G_3+\frac{s\cdot F_{3}^{2}}{G_4+\frac{s\cdot F_{4}^{2}}{G_5+\frac{\ddots }{\ddots G_n+\frac{s\cdot F_{n}^{2}}{G_{n+1}}}}}}}}\)

where \(G_n\) are the terms of a Fibonacci-like sequence, such that \(G_1=F_{n-1}\) and \(G_2=F_{n-1}+F_{n+2}\). That's how I've been doing it: the two last elements of both sequences are calculated and the sums are computed starting from the last terms of the series and of the continued fraction.

Now, why waiting 60 milliseconds when 30 millisecond will do for 1K digits? :-)



Same as the previous method, except that two terms of the series and two terms of the continued fraction are processed at a time.

OPTION ARITHMETIC DECIMAL_HIGH
! Reciprocal Fibonacci Constant
DEF SQRT(x) = EXP(LOG(x)/2)
DEF Fib(x) = INT((phi^x - (-phi)^(-x))/SQRT(5) + 0.01)
INPUT PROMPT "Number of digits: ":f
LET t = TIME
LET phi = (1 + SQRT(5))/2
LET k = CEIL(SQRT(f)*(phi - 1/(5*SQRT(SQRT(SQRT(f))))) - 1/2)
LET k = k + MOD(k,2) ! if k is odd then k = k + 1
LET sr = 0
LET a = Fib(k)
LET b = Fib(k - 1)
LET cf = 0
LET d = (Fib(k - 1) + Fib(k + 2))*Fib(k) + Fib(k - 1)^2
LET e = Fib(k - 1)*(Fib(k - 2) + Fib(k - 1) + Fib(k + 2))
FOR i = 1 TO k/2
LET sr = sr + (a + b)/(a*b) ! or sr = sr + 1/a + 1/b
LET cf = b*b*(cf + d)/(a*a - e*(cf + d)) ! faster alternative to
LET a = a - b ! cf = -b*b/(e - a*a/(cf + d))
LET b = b - a
LET d = d - e
LET e = e - d
NEXT i
LET cf = 1/(cf + d)
LET r$ = STR$(sr + cf)
LET r = TIME - t
PRINT
PRINT r$(1:2);
FOR i = 3 TO f + 1
PRINT r$(i:i);
IF MOD((i - 2),10) = 0 THEN PRINT " ";
IF MOD((i - 2),50) = 0 THEN
PRINT
PRINT " ";
END IF
NEXT i
IF MOD (i - 3,50) <> 0 OR f = 0 THEN PRINT
PRINT
PRINT "Runtime: ";
PRINT USING "0.##": r;
PRINT " seconds"
END

I may provide an RPL version later, but I will update this post rather than adding a new one.

-------------------------------------------------------------------------------------------

( 10-31-2018 12:59 PM ) Update:

RPL version

%%HP: T(3)A(R)F(.);
\<< PUSH RAD -105 CF -3 CF DUP \v/ DUP \v/ \v/ 5 * INV NEG 5 \v/ 1 + 2 / + * 2 INV - CEIL DUP 2 MOD + DUP 0 ROT
[[ 1 1 ]
 [ 1 0 ]] SWAP DUP2 ^ 3 GETI UNROT GET 4 ROLLD 4 ROLLD DUP + ^ 1 GETI UNROT GET 1 + 0 1 8 ROLL 2 /
  START PICK3 + DUP PICK3 * NEG 6 PICK SQ + / 4 PICK SQ * EXPAND ROT PICK3 - ROT OVER - ROT 6 ROLL 6 ROLL 6 ROLL + LASTARG * LASTARG 5 ROLLD 5 ROLLD / + ROT PICK3 - ROT OVER - 6 ROLL 6 ROLL 6 ROLL
  NEXT ROT + INV 5 ROLL + EXPAND 4 ROLLD 3 DROPN FXND PICK3 ALOG OVER - PICK3 * SWAP IQUOT + \->STR DUP HEAD -51 FC? { "." } { "," } IFTE + SWAP TAIL + 1 ROT 2 + SUB POP
\>>

Checksum #FF1Dh, 464.5 bytes

1000 decimal digits in 22m23.4s on the HP 50g

5000 decimal digits in 22m03.1s on the emulator (@ 2.4GHz)

3.
35988566624317755317201130291892717968890513373196848649555381532513031899668338​36154162164567900872
97045342928853913304136789017100883679591351733077119078580333550332507753187599​85048717977789700603
95645092153758927752656733540240331694417992939346109926262579646476518686594497​10216558984360881472
69324959107947387367337852332687749976272775794685367691854198146766874299876738​20969139012177220244
05208151094264934951374541667278955344470777775847802596340769074847415557910420​06750152034107053352
85129792635242062267537568055761955669720848843854407983324292851368070827522662​57975118864646409673
74615723872362955620536122030246354092526784242243470363103632014662980402490155​78724456176000319551
98790596994202917886694917480809674652368265408693839906987321175216695706385941​18145536473642687824
62926166650100098903804823359519893146150108288726392887669917149304053057745574​32156116729898561772
97313953707352919668843278980221650475850280918062910024442770174602410404177860​69190065037142835294
24549060159304435523683923718504987198630763391594936596142988354724590992947082​67223120172721663872
12418231942461961057752016152749577881832442620171321193355424014211322226700585​07790336805573465162
50584022990863875041525746936821124718256844645556016316277134557143751455007624​87099156182826892640
77299815147332688964961557144922458658613241975424930678178625935855419053397178​91505455586364455620
79380314058450581455450721741553866439013072379329369332459340804467686577840688​80825674051153703481
65653540656305697845332620804990536834597514423561119839926477047427769803101105​88155817184923932868
94389844332490569195491726876095595080412484819283259716387308029253819599769804​15448551217281602262
03319430666297327595319222766597355323536048195252102867001421399797884243553601​67959912204261837046
15748040151050250099917594368466922195960627316577648767941480521524363689518527​44138140626600534344
11601309785072730856410374536499056449092564022479235893154871766949243670359504​30223730003069311959
15230842412942958779782650009370949482395881794956517579528595425793938129592869​01360210094380436689
52309699341361290441703462826417037839835745622337038871970097160076557065423903​88747191135079115407
61134605405981579190230575626467157535916491877701487554582185026269324924055408​52156225044719698236
66233554986281689696668362778907832244552306778661186190701466995901904809618163​98991488427181733025
77991640047704791863450726856393944964259507568329203695785873826090877061694844​35450991486528604325
65638390016063340306115630353963540983632188358344730808839012965862794845736191​92828722289336977367
97823787549779581928084215139848092594370773749900079084904560712288744500834339​34927866138800198957
34867168389420674622442797718116543703036057676796294039383792686684190008693733​39749628021473634107
07806920809164559601268404152602112999654043604649332695667210328259257153518130​46845293343263071155
37801746317153246120127974550540292940747880846608451924758994639761595127788430​13427331220165756756
58061486335392946951764418852180537667074565698606044208952610950854918737031740​92680684538011907540
12913299974871973993118192889693921813655838683517341556816784000346983197133460​50391102634546299985
65647572232560931118513656555744824808892274072469603237445410557855544407201664​67493389833937288915
12205064177914205130750471784644944811905803552338083847929570164900961349168885​81297171817768451298
16652270445878950643816324667166234569860075170704208767090542338387590706535730​07253574394206705393
51489016774684898023581685604044328368996251110689026199698654113958550047889172​41289286530573960459
52956304818984785977933152007041340701120950340027986067411844908423122100242875​25516602104585306029
22892238520061424735392250410718191969280568460637845329235910300254887034488027​27717215706091545941
74746438950439128108937317098157904566980599145151178222281338283899824235972028​27748945586513936391
34524838564977794539084165164609787529500174427547604478851856024688447345148360​66921399141711901083
16900711733876744183059049274098365820972649176175217705801716875998786866072791​50571455108080357714
74910503685344257889044201305572642593686973313944115071648015446749503585463776​07632431556042666304
63955425809032409135700905265248048300624653066785950832462808429589748353027890​09884602603688205927
58052463208772411758069195877758240303454715751872910489968583866531493380284865​36933731835469634484
15044092149892251699966494897284989778964037755628480493187519547499568304516500​35737932722881001804
86457623483635176237516590209507046974062033317716447707406614289297307203270308​07528292464987612187
12329466370898509701781073945698974121734077615085286740959807507181162419388588​64474274231665018464
71554569249957276937832824696006445299378419273590386568422588998302821052023840​15914721662892536042
50632523916623289999832611104216389122078631278048053253887207877887998723652440​78850903523599082782
84141818003993773910950712076054580652250831335791926747402917363052309941864828​05359949788319437063


-------------------------------------------------------------------------------------------

Based primarily on the following DECIMAL BASIC version:

OPTION ARITHMETIC DECIMAL_HIGH
! Reciprocal Fibonacci Constant
DEF SQRT(x) =  EXP(LOG(x)/2)                                    
DEF Fib(x) = CEIL((phi^x - 1)/SQRT(5))
INPUT  PROMPT "Number of digits: ":f                           ! works for f <= 510 digits                                                
LET t = TIME                          
LET phi = (1 + SQRT(5))/2
LET k = CEIL(SQRT(f)*(phi - 1/(5*SQRT(SQRT(SQRT(f))))) - 1/2)
LET k = k + MOD(k,2)                                           ! if k is odd then k = k + 1
LET sr = 0                                                     ! (make sure k is even)
LET a = Fib(k) 
LET b = Fib(k - 1)
LET cf = 0
LET d = Fib(2*k + 1)                                           ! correct results for k <= 32 
LET e = Fib(2*k) + 1                                           ! for odd k, e = Fib(2*k) - 1 
FOR i = 1 TO k/2                                                 
   LET sr = sr + (a + b)/(a*b)                                 ! or sr = sr + 1/a + 1/b                  
   LET cf = b*b*(cf + d)/(a*a - e*(cf + d))                    ! faster alternative to 
   LET a = a - b                                               ! cf = -b*b/(e - a*a/(cf + d))
   LET b = b - a
   LET d = d - e
   LET e = e - d
NEXT i
LET cf = 1/(cf + d)
LET r$ = STR$(sr + cf)
LET r = TIME - t                            
PRINT
PRINT r$(1:2);
FOR i = 3 TO f + 1
   PRINT r$(i:i);
   IF MOD((i - 2),10) = 0 THEN PRINT " ";
   IF MOD((i - 2),50) = 0 THEN 
      PRINT
      PRINT "  ";
   END IF
NEXT i
IF MOD (i - 3,50) <> 0  OR f = 0 THEN PRINT
PRINT 
PRINT "Runtime: ";
PRINT  USING "0.##": r;
PRINT " seconds"
END

-------------------------------------------------------------------------------------------

Update:

Rewriting of formulae and exact expressions for n (number of terms) as function of d (number of correct significant digits).

\(\psi \approx \frac{1}{F_{1}}+\frac{1}{F_{2}}+\frac{1}{F_{3}}+\cdots +\frac{1}{F_{n-1}}+\frac{1}{F_{n}}+\frac{1}{F_{n-1}+\frac{(-1)^{n+1} F_{1}^{2}}{F_{n+2}\cdot F_{1}+F_{n-1}\cdot F_{2}+\frac{(-1)^{n+1} F_{2}^{2}}{F_{n+2}\cdot F_{2}+F_{n-1}\cdot F_{3}+\frac{(-1)^{n+1} F_{3}^{2}}{F_{n+2}\cdot F_{3}+F_{n-1}\cdot F_{4}+\frac{(-1)^{n+1} F_{4}^{2}}{F_{n+2}\cdot F_{4}+F_{n-1}\cdot F_{5}+\frac{\ddots }{\ddots F_{n+2}\cdot F_{n-1}+F_{n-1}\cdot F_{n}+\frac{(-1)^{n+1} F_{n}^{2}}{F_{n+2}\cdot F_{n}+F_{n-1}\cdot F_{n+1} }}}}}}}\)

or, more compactly,

\(\psi \approx \frac{1}{F_{1}}+\frac{1}{F_{2}}+\frac{1}{F_{3}}+\cdots +\frac{1}{F_{n-1}}+\frac{1}{F_{n}}+\frac{1}{G_1+\frac{(-1)^{n+1} F_{1}^{2}}{G_2+\frac{(-1)^{n+1} F_{2}^{2}}{G_3+\frac{(-1)^{n+1} F_{3}^{2}}{G_4+\frac{(-1)^{n+1} F_{4}^{2}}{G_5+\frac{\ddots }{\ddots G_n+\frac{(-1)^{n+1} F_{n}^{2}}{G_{n+1}}}}}}}}\)


where \(G_n\) are the terms of a Fibonacci-like sequence, such that \(G_1=F_{n-1}\) and \(G_2=F_{n-1}+F_{n+2}\).


The number of terms of the series plus the number of the terms of the continued fraction required to yield d correct significant digits is given by

\(n = \left \lceil \frac{1}{2}\sqrt{\frac{d\cdot \ln \left ( 100 \right ) + \ln \left ( 5 \right )}{\ln \left ( \varphi \right )}}-1 \right \rceil\)

where \(\varphi = \frac{1+\sqrt{5}}{2}\) is the golden ratio.


Other equivalent expressions are

\(n =\left \lceil \frac{1}{2}\sqrt{\frac{d\cdot \ln \left ( 100 \right ) + \ln \left ( 5 \right )}{ \sinh^{-1}\left ( \frac{1}{2} \right ) }}-1 \right \rceil\)

and

\(n =\left \lceil \frac{1}{2}\sqrt{\frac{2\left ( d-1 \right ) + \log_{10} \left ( 500 \right )}{\log_{10} \left ( \varphi \right ) }}-1 \right \rceil\)

RPL (HP-49G/G+/50g), 460.5 bytes, Checksum # 5904h

%%HP: T(3)A(R)F(.);
\<< RCLF SWAP RAD -105 CF -3 CF R\->I DUP 1 + 100 LN *
5 LN + 2 INV ASINH / \v/ 2 / CEIL DUP 2 MOD + DUP 0 ROT
[[ 1 1 ]
[ 1 0 ]] SWAP DUP2 ^ 3 GETI UNROT GET 4 ROLLD 4 ROLLD
DUP + ^ 1 GETI UNROT GET 1 + 0 1 8 ROLL 2 /
 START PICK3 + DUP PICK3 * NEG 6 PICK SQ + / 4 PICK SQ
* EXPAND ROT PICK3 - ROT OVER - ROT 6 ROLL 6 ROLL 6
ROLL + LASTARG * LASTARG 5 ROLLD 5 ROLLD / + ROT PICK3
- ROT OVER - 6 ROLL 6 ROLL 6 ROLL
 NEXT ROT + INV 5 ROLL + EXPAND 4 ROLLD 3 DROPN FXND
PICK3 ALOG OVER - PICK3 * SWAP IQUOT + \->STR DUP HEAD
-51 FC? { "." } { "," } IFTE + SWAP TAIL + 1 ROT 2 +
SUB SWAP STOF
\>>

DECIMAL BASIC

OPTION ARITHMETIC DECIMAL_HIGH
! Reciprocal Fibonacci Constant
DEF Fib(x) = CEIL((phi^x - 1)/SQR(5))
LET phi = (1 + SQR(5))/2
INPUT PROMPT "Number of decimal places: ":n
IF n < 500 THEN LET l = 0 ELSE LET l = 1
LET t = TIME                                                     ! k: # of terms of the series & cont'ed fraction 
LET k = CEIL(EXP(LOG((2*n + LOG10(500))/LOG10(phi))/2)/2 - l)    ! k = Sqrt((2*n + log10(500))/log10(phi))/2 - 1                                                                                                                                        
LET k = k + MOD(k,2)                                             ! if k is odd then k = k + 1
LET sr = 0                                                       ! (make sure k is even)
LET a = Fib(k) 
LET b = Fib(k - 1)
LET cf = 0
LET d = Fib(2*k + 1)                                            
LET e = Fib(2*k) + 1                                             ! for even k, e = Fib(2*k) + 1 
FOR i = 1 TO k/2                                                 !(for odd k, e = Fib(2*k) - 1)                                     
   LET sr = sr + (a + b)/(a*b)                                   ! or sr = sr + 1/a + 1/b                  
   LET cf = b*b*(cf + d)/(a*a - e*(cf + d))                      ! faster alternative to 
   LET a = a - b                                                 ! cf = -b*b/(e - a*a/(cf + d))
   LET b = b - a
   LET d = d - e
   LET e = e - d
NEXT i
LET cf = 1/(cf + d)
LET r$ = STR$(sr + cf)
LET r = TIME - t                            
PRINT
PRINT r$(1:2);
FOR i = 3 TO n + 2
   PRINT r$(i:i);
   IF MOD((i - 2),10) = 0 THEN PRINT " ";
   IF MOD((i - 2),50) = 0 THEN 
      PRINT
      PRINT "  ";
   END IF
NEXT i
IF MOD (i - 3,50) <> 0  OR n = 0 THEN PRINT
PRINT 
PRINT "Runtime: ";
PRINT  USING "0.##": r;
PRINT " seconds"
END


Number of decimal places: 1000

3.3598856662 4317755317 2011302918 9271796889 0513373196 
  8486495553 8153251303 1899668338 3615416216 4567900872 
  9704534292 8853913304 1367890171 0088367959 1351733077 
  1190785803 3355033250 7753187599 8504871797 7789700603 
  9564509215 3758927752 6567335402 4033169441 7992939346 
  1099262625 7964647651 8686594497 1021655898 4360881472 
  6932495910 7947387367 3378523326 8774997627 2775794685 
  3676918541 9814676687 4299876738 2096913901 2177220244 
  0520815109 4264934951 3745416672 7895534447 0777775847 
  8025963407 6907484741 5557910420 0675015203 4107053352 
  8512979263 5242062267 5375680557 6195566972 0848843854 
  4079833242 9285136807 0827522662 5797511886 4646409673 
  7461572387 2362955620 5361220302 4635409252 6784242243 
  4703631036 3201466298 0402490155 7872445617 6000319551 
  9879059699 4202917886 6949174808 0967465236 8265408693 
  8399069873 2117521669 5706385941 1814553647 3642687824 
  6292616665 0100098903 8048233595 1989314615 0108288726 
  3928876699 1714930405 3057745574 3215611672 9898561772 
  9731395370 7352919668 8432789802 2165047585 0280918062 
  9100244427 7017460241 0404177860 6919006503 7142835296 
  
Runtime: 0.01 seconds