Reciprocal Fibonacci constant

Question

Seeing as there have been an awful lot of normal Fibonacci challenges, I decided that it might be interesting to calculate the Reciprocal Fibonacci constant - namely, the sum of the reciprocals of the Fibonacci sequence.

The challenge is to calculate the Reciprocal Fibonacci constant with the number of Fibonacci series to use digit given as input, i.e. an input of 10 means to calculate based on the reciprocals of the first 10 Fibonacci numbers. In the likely case of a tie, shortest code wins.

The winner will be chosen in three weeks via standard code-golf rules.

Answer 1

K (19)

(or 17 if you skip defining it as a function)

f:{+/*+%x(|+\)\|!2}

Tested on Kona.

Basically, it generates the first x values of the fibonacci sequence into an array (without using builtins), divides 1 by each value of the entire array, transposes and sums it up.

(thanks to @tmartin for the better sum method)

Answer 2

Perl - 35 bytes

print!map$\-=1/($%+=$.=$%-$.),0..<>

Sample usage:

$ echo 10 | perl inv-fib-sum.pl
3.34170499581934

---

Further Analysis

It's interesting to note that the sum

is convergent. Supposing we wanted to calculate a few thousand digits or so, the naïve approach is almost sufficient. The convergence is quite slow at first, but speeds up rapidly, so that 1000 digits only takes about 4800 terms. A sample Python implementation might be:

a=[1,1]
for i in range(4800):a=[a[0]+a[1]]+a
z=10**1000
print sum(map(lambda i:z/i,a))

which after a second or so outputs:

33598856662431775531720113029189271796889051337319684864955538153251303189966833836154162164567900872970453429288539133041367890171008836795913517330771190785803335503325077531875998504871797778970060395645092153758927752656733540240331694417992939346109926262579646476518686594497102165589843608814726932495910794738736733785233268774997627277579468536769185419814676687429987673820969139012177220244052081510942649349513745416672789553444707777758478025963407690748474155579104200675015203410705335285129792635242062267537568055761955669720848843854407983324292851368070827522662579751188646464096737461572387236295562053612203024635409252678424224347036310363201466298040249015578724456176000319551987905969942029178866949174808096746523682654086938399069873211752166957063859411814553647364268782462926166650100098903804823359519893146150108288726392887669917149304053057745574321561167298985617729731395370735291966884327898022165047585028091806291002444277017460241040417786069190065037142832933

(The last four digits don't quite converge, but we'll ignore that for now.)

Let's try to speed up the convergence a bit. A standard trick is to use Euler's Transform. After expansion and simplification, this produces a nicer result:

It should be fairly apparent why this converges more quickly; each term has 3 terms in the denominator rather than just one. Calculating 1000 digits takes only 1600 (one third as many) terms:

a=[1,1]
for i in range(1601):a=[a[0]+a[1]]+a
z=10**1000
print sum(map(lambda i:(-1)**i*z/(a[i]*a[i+1]*a[i+2]),range(1601)))

Output:

3598856662431775531720113029189271796889051337319684864955538153251303189966833836154162164567900872970453429288539133041367890171008836795913517330771190785803335503325077531875998504871797778970060395645092153758927752656733540240331694417992939346109926262579646476518686594497102165589843608814726932495910794738736733785233268774997627277579468536769185419814676687429987673820969139012177220244052081510942649349513745416672789553444707777758478025963407690748474155579104200675015203410705335285129792635242062267537568055761955669720848843854407983324292851368070827522662579751188646464096737461572387236295562053612203024635409252678424224347036310363201466298040249015578724456176000319551987905969942029178866949174808096746523682654086938399069873211752166957063859411814553647364268782462926166650100098903804823359519893146150108288726392887669917149304053057745574321561167298985617729731395370735291966884327898022165047585028091806291002444277017460241040417786069190065037142834500

(Here again, the last 4 digits don't converge.)

We're not quite done yet. If we combine adjacent terms, we end up with the following:

Factoring out each term from the remainder of the summation gives the nested expression:

Now we're getting somewhere. Notice that the numerators of follow OEIS A206351 (with the exception of the first term, which is doubled):

and the denominators follow OEIS A081016 (shifted by one term):

Each of these have very simple recurrence relations, namely:

and

respectively. Putting it all together, we find that we need only 800 iterations for 1000 digits:

b,c=[16,3,1],[273,40,3]
for i in range(800):b,c=[7*b[0]-b[1]-4]+b,[7*c[0]-c[1]-1]+c
s=z=10**1000
for x,y in zip(b,c):s=(z+s)*x/y
print s

which outputs:

3598856662431775531720113029189271796889051337319684864955538153251303189966833836154162164567900872970453429288539133041367890171008836795913517330771190785803335503325077531875998504871797778970060395645092153758927752656733540240331694417992939346109926262579646476518686594497102165589843608814726932495910794738736733785233268774997627277579468536769185419814676687429987673820969139012177220244052081510942649349513745416672789553444707777758478025963407690748474155579104200675015203410705335285129792635242062267537568055761955669720848843854407983324292851368070827522662579751188646464096737461572387236295562053612203024635409252678424224347036310363201466298040249015578724456176000319551987905969942029178866949174808096746523682654086938399069873211752166957063859411814553647364268782462926166650100098903804823359519893146150108288726392887669917149304053057745574321561167298985617729731395370735291966884327898022165047585028091806291002444277017460241040417786069190065037142835294

(Here, finally, the last 4 digits converge correctly.)

But that's still not quite everything. If we observe the intermediate values for s, we find that it converges to a different value entirely before converging on the actual convergence point. The reason for this is the following:

Solving for a stable s, we find that:

Because this is a simple root, we can use Newton's Method to get us most of the way there, and then jump in at a much later point in the iteration. Only about 400 digits of precision are necessary (as the b and c values aren't any larger than that anyway), which can be achieved with just 7 iterations, while saving 320 iterations of the main loop:

b,c=[16,3,1],[273,40,3]
for i in range(480):b,c=[7*b[0]-b[1]-4]+b,[7*c[0]-c[1]-1]+c
z=10**1000;s=z/17
for i in range(7):s-=(s*s+s*z-z*z/16)/(2*s+z)
for x,y in zip(b,c):s=(z+s)*x/y
print s

Output is identical to the previous, runtime is about 0.02s on my system using PyPy v2.1. Even though it needs one tenth the number of iterations as the original, it's significantly faster than 10x due to multiplying and dividing by much smaller terms. I don't think much more can be tweaked out of it, although I'd be happy to be shown wrong.

Answer 3

Mathematica (32 characters without built-in Fibonacci)

Tr[1/⌈(.5+√5/2)^Range@#/√5-.5⌉]&

Here I used the rounding formula for Fibonacci numbers

where φ is the golden ratio.

Answer 4

Kona (25 21)

{+/%(x(|+\)\1 1)[;1]}

Probably could be made smaller by experts, but I am still learning the language.

  f 10
3.341705
  f 3
2.8333
  f 25
3.359872
  f 400
3.359886

The last one didn't actually take any more time than the others.

Answer 5

Mathematica 30 24

With 6 characters saved thanks to ybeltukov.

 Tr[1/Fibonacci@Range@#]&

---

Before adding the terms:

1/Fibonacci@Range@#&[20]

---

With addition included:

 Tr[1/Fibonacci@Range@#]&[20]

Answer 6

APL, 21 chars/bytes*

{+/÷{↑+\∘⌽⍣⍵⊢0 1}¨⍳⍵}

Example

      {+/÷{↑+\∘⌽⍣⍵⊢0 1}¨⍳⍵} 10
3.330469041
      ⍪{+/÷{↑+\∘⌽⍣⍵⊢0 1}¨⍳⍵}¨⍳10
1          
2          
2.5        
2.833333333
3.033333333
3.158333333
3.23525641 
3.282875458
3.312287223
3.330469041
      {+/÷{↑+\∘⌽⍣⍵⊢0 1}¨⍳⍵} 100
3.359885666

⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯
^{*: APL can be written in its own (legacy) single-byte charset that maps APL symbols to the upper 128 byte values. Therefore, for the purpose of scoring, a program of N chars that only uses ASCII characters and APL symbols can be considered to be N bytes long.}

Answer 7

Javascript, 33

Input: n = 10

s=f=t=1;for(;n--;)t+=1/(s+=f=s-f)

Based in my first post in Perl, here is the result directly from the Google Chrome Console.

Answer 8

K, 22

...

{+/%x#x{x,+/-2#x}/1 1}

Answer 9

Perl 6, 27 bytes

{sum 1 X/(1,&[+]...*)[^$_]}

Try it online!

Alternatively

{sum 1 xx$_ Z/(1,&[+]...*)}

Try it online!

Answer 10

GTB, 35

Runs on a TI-84 calculator

1→Y:0→X`N4;N,1,N~1/X:Y→Z:X+Y→Y:Z→X&

• Store 1 in Y and 0 in X
• Take N as input
• Initialize a For loop
• Display X
• Store Y in Z and X+Y in Y
• Store Z in X
• End the For loop

Answer 11

PERL, 62 43

$s=$t=1;$t+=1/($a=$f+$s),$f=$s,$s=$a,for 0..<>;say $t

Edit:

After some more time looking at my code I was able to save 19 chars:

$s=1;$t+=1/($s+=$f=$s-$f)for 0..<>;say $t+2

input: 10
> 3.35294128575734

Answer 12

BC - 116

define f(n){if(n < 2){return n;}else{return f(n-1)+f(n-2);}}
define r(n){for(i=1;i<n;++i){m=m+(1/f(i));};return m;}

call r(n) with n to solve. The precision can be set arbitrarily, so this is the most accurate solution.

Answer 13

Forth, 64

: r 1 1 rot 1e 0 do 2dup - nip tuck + dup s>f 1/f f+ swap loop ;

usage:

10 r f. 3.34170499581934  ok
100 r f. 3.35988566624318  ok
1000 r f. 3.35988566624318  ok

Answer 14

Python (55 51)

i,r=p=1,1
while i<n+1:r+=1./p[-1];p=p[1],sum(p);i+=1
r

i,r,a,b=[1]*4
while i<n+1:r+=1./b;a,b=b,b+a;i+=1
r

In interactive mode:

n=10
3.341704995819338

n=100

3.3598856429940778

3.359885666243178

n=400

3.3598855578800064

3.359885666243178

Answer 15

RPL (HP 49G/G+/50g), 50 bytes

« 0. 1. DUP2 5. ROLL
  START SWAP OVER + ROT OVER INV + UNROT
  NEXT DROP2
»

Examples:

10 -> 3.33046904076

11 -> 3.34170499582

30 -> 3.35988372158

55 -> 3.35988566622

In theory, the sum of the first 55 digits would give 12 correct digits. The last digit should be '4' instead of '2'. This should be fixed by adding up the terms backwards, but the code would be slightly longer.

If the constant were to be computed to 1000 decimal places using the definition the sum should be carried out to at least the first 4790 terms, which would take too much time on the HP 50g, if possible at all. For this task a more efficient method is required, like the one used in the following RPL program (464 bytes, checksum # 207Eh) whose argument is the desired number of decimal places:

%%HP: T(3)A(R)F(.);
\<< PUSH RAD -105 CF -3 CF R\->I DUP 1 + 2 INV SWAP 100 LN * 5 LN + OVER ASINH / \v/ * 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
\>>

'RFC' STO

1000 RFC -->

3.3598856662431775531720113029189271796889051337319684864955538153251303189966833836154162164567900872970453429288539133041367890171008836795913517330771190785803335503325077531875998504871797778970060395645092153758927752656733540240331694417992939346109926262579646476518686594497102165589843608814726932495910794738736733785233268774997627277579468536769185419814676687429987673820969139012177220244052081510942649349513745416672789553444707777758478025963407690748474155579104200675015203410705335285129792635242062267537568055761955669720848843854407983324292851368070827522662579751188646464096737461572387236295562053612203024635409252678424224347036310363201466298040249015578724456176000319551987905969942029178866949174808096746523682654086938399069873211752166957063859411814553647364268782462926166650100098903804823359519893146150108288726392887669917149304053057745574321561167298985617729731395370735291966884327898022165047585028091806291002444277017460241040417786069190065037142835294

(22 min 23 sec, on the real HP 50g)

---

For one thousand digits the first 50 terms of the series plus another 50 terms of a continued fraction are evaluated, two at a time in only 25 iterations (48 terms of each would have sufficed, though). An equivalent DECIMAL BASIC program takes just 10 milliseconds on my computer (Intel(R) Core(TM) Duo CPU E4700 @ 2.6 GHz).

Answer 16

Japt -mx, 6 bytes

1/MgUÄ

Try it

Answer 17

JAEL, 9 7 bytes

Changing the diacritics from one character to another gave me 1 byte

Still working on SBCS encoding.

#`&@uȦ

Explanation (generated automatically):

./jael -u explain '#`&@uȦ'

ORIGINAL CODE:
#`&@uȦ

EXPANDING EXPLANATION:
Ȧ => .a!

EXPANDED CODE:
#`&@u.a!

COMPLETED CODE:
#`&@u.a!,


#       ,               repeat (p1) times:
 `                              stack 1
  &                             push number of iterations of this loop
   @                            push nth fibonacci number
    u                           push p2 / p1
     .                          push the value under the tape head
      a                         push p1 + p2
       !                        write p1 to the tapehead
         ␄              print machine state

Answer 18

C# (.NET Core), 66 bytes

a=>{double r=1,x=1,y=1,i=0;for(;++i<a;y+=x,x=y-x)r+=1/y;return r;}

Try it online!

Can't use floats because of imprecision.

Example where input = 8:

• double: 3.28287545787546 (rounds to 3.282875)
• float: 3.282876 (off by 0.000001)

Ungolfed:

a => {
    double r = 1,                       // initialize r (sum of reciprocals) - start at 1 to save one byte in for loop
            x = 1, y = 1,                   // initialize x (2nd largest Fibonacci) and y (largest Fibonacci)
            i = 0;                          // initialize i (loop counter)
    for(; ++i < a; y += x, x = y - x)   // from 1 to a, while incrementing the Fibonacci numbers
        r += 1 / y;                         // increment the reciprocal Fibonacci
    return r;                           // return the reciprocal Fibonacci
}

Answer 19

cQuents, 13 11 bytes

;/b$
=1:Z+Y

Try it online!

Explanation

;          The sum of the first n terms of the sequence
 /         which are 1 divided by
  b$       the nth term in the sequence on the second line

=1:Z+Y     Fibonacci with starting indexes 1,1

Answer 20

Jelly, 5 bytes

RÆḞİS

Try it online!

How it works

RÆḞİS    Main link (monad). Input: integer
R        Range
 ÆḞ      nth Fibonacci number
   İ     Reciprocal
    S    Sum