Search
Member List
Calendar
Help
|
Mardi Gras True Fibs
|
|
03-01-2017, 12:41 AM
(This post was last modified: 03-01-2017 01:55 AM by Gerson W. Barbosa.)
Post: #1
|
|||
|
|||
|
Mardi Gras True Fibs
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: 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.359885666243177553172011302918927179688905133731968486495553815325130318996683383615416216456790087 2 12 3.359885666242129988464989239268857623668736813180767663357400070297267494293999199307997271121544920 2 20 3.359885666243177553167433281787748928648776577928290964879706040408646726702491028523646222602832650 2 30 3.359885666243177553172011302918764514065619978389958675293504780890783041151282869607821301725790765 2 40 3.359885666243177553172011302918927179688899353919213085077297503281892745245057092283592271334021082 2 50 3.359885666243177553172011302918927179688905133731968281128034773312024205107370483101951938492442443 4 10 3.359885666243177552933037307174019435890442521909028053549830953375137979954347559962927929423499662 4 20 3.359885666243177553172011302918927179651795581225695045861072258099578065882672826328709609425804709 4 30 3.359885666243177553172011302918927179688905133731968486489791485719247471989816805452564095340121083 4 40 3.359885666243177553172011302918927179688905133731968486495553815325130318995788615553566782514358382 4 50 3.359885666243177553172011302918927179688905133731968486495553815325130318996683383615416216456651148 10 20 3.359885666243177553172011302918927179688905133731968486495553815325130318996683383615416201795067149 11 20 3.359885666243177553172011302918927179688905133731968486495553815325130318996683383615416216456790047 14 14 3.359885666243177553172011302918927179688905133731968486495553815325130318996683383615414413089430501 15 14 3.359885666243177553172011302918927179688905133731968486495553815325130318996683383615416216455221736 14 16 3.359885666243177553172011302918927179688905133731968486495553815325130318996683383615416216456790085 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: 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. |
|||
|
|
03-01-2017, 01:32 AM
Post: #2
|
|||
|
|||
|
RE: Mardi Gras True Fibs
Gerson, you continue to amaze!
|
|||
|
|
|
03-01-2017, 04:34 AM
(This post was last modified: 03-01-2017 05:29 AM by Gerson W. Barbosa.)
Post: #3
|
|||
|
|||
|
RE: Mardi Gras True Fibs
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.359885666243177553172011302918927179688905133731968486495553815325130318996683383615416216456790087297045342928853913304136789017100883679591351733077119078580333550332507753187599850487179777897006039564509215375892775265673354024033169441799293934610992626257964647651868659449710216558984360881472693249591079473873673378523326877499762727757946853676918541981467668742998767382096913901217722024405208151094264934951374541667278955344470777775847802596340769074847415557910420067501520341070533528512979263524206226753756805576195566972084884385440798332429285136807082752266257975118864646409673746157238723629556205361220302463540925267842422434703631036320146629804024901557872445617600031955198790596994202917886694917480809674652368265408693839906987321175216695706385941181455364736426878246292616665010009890380482335951989314615010828872639288766991714930405305774557432156116729898561772973139537073529196688432789802216504758502809180629100244427701746024104041778606919006503714283529 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. |
|||
|
|
|
03-04-2017, 03:34 PM
(This post was last modified: 12-03-2017 09:49 AM by Gerson W. Barbosa.)
Post: #4
|
|||
|
|||
|
RE: Mardi Gras True Fibs
Here is a UserRPL solution that works on the HP 49G/50g and doesn't rely on external arbitrary precision libraries:
Code: 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.359885666243177553172011302918927179688905133731968486495553815325130318996683383615416216456790087297045342928 (112 digits, 13.40 s) 50 d --> 3.3598856662431775531720113029189271796889051337319 (6.00 s) 100 d --> 3.359885666243177553172011302918927179688905133731968486495553815325130318996683383615416216456790087 (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. |
|||
|
|
|
06-14-2018, 12:41 AM
(This post was last modified: 06-14-2018 12:55 AM by Gerson W. Barbosa.)
Post: #5
|
|||
|
|||
|
RE: Mardi Gras True Fibs
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 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 " "; END IF NEXT i IF MOD (i - 3,50) <> 0 OR f = 0 THEN 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 |
|||
|
|
|
06-17-2018, 06:00 PM
(This post was last modified: 11-12-2018 01:43 AM by Gerson W. Barbosa.)
Post: #6
|
|||
|
|||
RE: Mardi Gras True Fibs
(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} }}}}}}}\) 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? :-) ![[Image: 41047424120_d86a38eb1e_b.jpg]](https://c2.staticflickr.com/2/1746/41047424120_d86a38eb1e_b.jpg) 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 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 " "; END IF NEXT i IF MOD (i - 3,50) <> 0 OR f = 0 THEN 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. 3598856662431775531720113029189271796889051337319684864955538153251303189966833836154162164567900872 9704534292885391330413678901710088367959135173307711907858033355033250775318759985048717977789700603 9564509215375892775265673354024033169441799293934610992626257964647651868659449710216558984360881472 6932495910794738736733785233268774997627277579468536769185419814676687429987673820969139012177220244 0520815109426493495137454166727895534447077777584780259634076907484741555791042006750152034107053352 8512979263524206226753756805576195566972084884385440798332429285136807082752266257975118864646409673 7461572387236295562053612203024635409252678424224347036310363201466298040249015578724456176000319551 9879059699420291788669491748080967465236826540869383990698732117521669570638594118145536473642687824 6292616665010009890380482335951989314615010828872639288766991714930405305774557432156116729898561772 9731395370735291966884327898022165047585028091806291002444277017460241040417786069190065037142835294 2454906015930443552368392371850498719863076339159493659614298835472459099294708267223120172721663872 1241823194246196105775201615274957788183244262017132119335542401421132222670058507790336805573465162 5058402299086387504152574693682112471825684464555601631627713455714375145500762487099156182826892640 7729981514733268896496155714492245865861324197542493067817862593585541905339717891505455586364455620 7938031405845058145545072174155386643901307237932936933245934080446768657784068880825674051153703481 6565354065630569784533262080499053683459751442356111983992647704742776980310110588155817184923932868 9438984433249056919549172687609559508041248481928325971638730802925381959976980415448551217281602262 0331943066629732759531922276659735532353604819525210286700142139979788424355360167959912204261837046 1574804015105025009991759436846692219596062731657764876794148052152436368951852744138140626600534344 1160130978507273085641037453649905644909256402247923589315487176694924367035950430223730003069311959 1523084241294295877978265000937094948239588179495651757952859542579393812959286901360210094380436689 5230969934136129044170346282641703783983574562233703887197009716007655706542390388747191135079115407 6113460540598157919023057562646715753591649187770148755458218502626932492405540852156225044719698236 6623355498628168969666836277890783224455230677866118619070146699590190480961816398991488427181733025 7799164004770479186345072685639394496425950756832920369578587382609087706169484435450991486528604325 6563839001606334030611563035396354098363218835834473080883901296586279484573619192828722289336977367 9782378754977958192808421513984809259437077374990007908490456071228874450083433934927866138800198957 3486716838942067462244279771811654370303605767679629403938379268668419000869373339749628021473634107 0780692080916455960126840415260211299965404360464933269566721032825925715351813046845293343263071155 3780174631715324612012797455054029294074788084660845192475899463976159512778843013427331220165756756 5806148633539294695176441885218053766707456569860604420895261095085491873703174092680684538011907540 1291329997487197399311819288969392181365583868351734155681678400034698319713346050391102634546299985 6564757223256093111851365655574482480889227407246960323744541055785554440720166467493389833937288915 1220506417791420513075047178464494481190580355233808384792957016490096134916888581297171817768451298 1665227044587895064381632466716623456986007517070420876709054233838759070653573007253574394206705393 5148901677468489802358168560404432836899625111068902619969865411395855004788917241289286530573960459 5295630481898478597793315200704134070112095034002798606741184490842312210024287525516602104585306029 2289223852006142473539225041071819196928056846063784532923591030025488703448802727717215706091545941 7474643895043912810893731709815790456698059914515117822228133828389982423597202827748945586513936391 3452483856497779453908416516460978752950017442754760447885185602468844734514836066921399141711901083 1690071173387674418305904927409836582097264917617521770580171687599878686607279150571455108080357714 7491050368534425788904420130557264259368697331394411507164801544674950358546377607632431556042666304 6395542580903240913570090526524804830062465306678595083246280842958974835302789009884602603688205927 5805246320877241175806919587775824030345471575187291048996858386653149338028486536933731835469634484 1504409214989225169996649489728498977896403775562848049318751954749956830451650035737932722881001804 8645762348363517623751659020950704697406203331771644770740661428929730720327030807528292464987612187 1232946637089850970178107394569897412173407761508528674095980750718116241938858864474274231665018464 7155456924995727693783282469600644529937841927359038656842258899830282105202384015914721662892536042 5063252391662328999983261110421638912207863127804805325388720787788799872365244078850903523599082782 8414181800399377391095071207605458065225083133579192674740291736305230994186482805359949788319437063 ------------------------------------------------------------------------------------------- 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 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 " "; END IF NEXT i IF MOD (i - 3,50) <> 0 OR f = 0 THEN 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 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 " "; END IF NEXT i IF MOD (i - 3,50) <> 0 OR n = 0 THEN 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 |
|||
|
|
|
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 1 Guest(s)