Mardi Gras True Fibs

[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