# Reciprocal Fibonacci constant

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.

• This is equal (if I've understood it right) to 1/φ (reciprocal of golden ratio). If you want us to actually use the Fibonacci numbers in the calculation, you should specify. If not, there's certainly languages where φ is a builtin. (not APL for a change) – marinus Jan 9 '14 at 22:40
• @marinus Changed. – Hosch250 Jan 9 '14 at 22:55
• @marinus, it's not equal to 1/phi. It does have a closed form, but it's quite tricky. Mathematica probably has a built-in, but I doubt many other languages do. – Peter Taylor Jan 9 '14 at 23:22
• @OP, "highest accuracy possible" is a useless winning criterion. Anyone who has a language which supports arbitrary decimal precision, or who can be bothered to write the support for it, can make an implementation with a precision parameter and then engage in an edit war to increase that parameter. It would make more sense to ask for a function which takes the precision as a parameter. Judging on speed is also tricky, because it depends on many factors (CPU model, RAM available, system load, ...). – Peter Taylor Jan 9 '14 at 23:28
• @PeterTaylor Is this better? – Hosch250 Jan 9 '14 at 23:41

## 20 Answers

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)

• A slight variation on yours for 17 chars -- {+/*+%x(|+\)\|!2} – tmartin Jan 13 '14 at 11:48
• @tmartin Thanks, yes, that was a less convoluted sum method :) – Joachim Isaksson Jan 13 '14 at 12:02

## 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+a]+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+a]+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-b-4]+b,[7*c-c-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-b-4]+b,[7*c-c-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.

# 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.

## 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.

• I am learning some languages by completing challenges in them. Nice to see that I am not the only one. – Hosch250 Jan 10 '14 at 3:20
• Save yourself two characters by just calling the function as a lambda rather than a named function. Then save a further two by using % as the reciprocal function rather than doing 1.% -- {+/%(x(|+\)\1 1)[;1]} for 21 chars – tmartin Jan 13 '14 at 11:18
• @tmartin: Odd...I thought I tried doing reciprocal and was getting 0 because of integer division, but it is currently working for me. I'm updating the answer, thanks a lot! – Kyle Kanos Jan 13 '14 at 14:03
• Are you sure you're computing the right sum? I get 2.833333333 for n=4, for example, instead of 3. One of us is wrong! – Tobia Jan 13 '14 at 18:15
• @Tobia: A difference between 0 and 1 indexing. f 0 results in 1.0, f 1 -> 2.0 and so on. – Kyle Kanos Jan 13 '14 at 18:37

# Mathematica 30 24

With 6 characters saved thanks to ybeltukov.

 Tr[1/Fibonacci@Range@#]&


Before adding the terms:

1/Fibonacci@Range@#& With addition included:

 Tr[1/Fibonacci@Range@#]& • Fibonacci is Listable :) It can be reduced to Tr[1/Fibonacci@Range@#]& – ybeltukov Jan 10 '14 at 0:57

# 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.

# 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. # K, 22

...

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


## GTB, 35

Runs on a TI-84 calculator

1→Y:0→XN4;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
• I assume it works, but can you provide an explanation? – Hosch250 Jan 10 '14 at 0:09
• I've added one. – Timtech Jan 10 '14 at 0:10

# 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.

# 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

## 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


## Python (55 51)

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


i,r,a,b=*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


# Perl 6, 27 bytes

{sum 1 X/(1,&[+]...*)[^$_]}  Try it online! Alternatively {sum 1 xx$_ Z/(1,&[+]...*)}


Try it online!

# Japt -mx, 6 bytes

1/MgUÄ


Try it

# 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
}


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

# 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

• Yes, you're right. I miscounted it. – Eduardo Hoefel Nov 8 '18 at 11:23
• I know the language is optimised for characters, but this site scores by bytes, so it's misleading to include the character count – Jo King Nov 8 '18 at 12:04
• You can get 8 bytes if you don't compress u.. On the topic of an SBCS, you might have trouble because your docs list over 600 commands or combinations of commands, and a SBCS language can only have 256 single byte commands. On the other hand, there's a large number of useless one to one translations... – Jo King Nov 8 '18 at 12:23
• Yes @JoKing, you are right. I started the whole idea a while ago and didn't consider that it would affect the byte size. Now I'm writing my undergraduate thesis on top of JAEL and it's too late to reconsider some decisions (I'm presenting it by the end of november). About the command list: yes, it was automatically generated and I didn't give it enough attention. As this is my first project on this area, I'll gather experience and try to learn from the bad decisions. Thank you for sharing your thoughts! – Eduardo Hoefel Nov 8 '18 at 12:54

# 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


# 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
`