Continued fraction: Difference between revisions

Content added Content deleted

Inline

Revision as of 09:30, 6 September 2012

A number may be represented as a continued fraction (see Mathworld for more information) as follows:

a_{0}+{\cfrac {b_{1}}{a_{1}+{\cfrac {b_{2}}{a_{2}+{\cfrac {b_{3}}{a_{3}+\ddots }}}}}}

The task is to write a program which generates such a number and prints a real representation of it. The code should be tested by calculating and printing the square root of 2, Napier's Constant, and Pi, using the following coefficients:

For the square root of 2, use $a_{0}=1$ then $a_{N}=2$ . $b_{N}$ is always $1$ .

{\sqrt {2}}=1+{\cfrac {1}{2+{\cfrac {1}{2+{\cfrac {1}{2+\ddots }}}}}}

For Napier's Constant, use $a_{0}=2$ , then $a_{N}=N$ . $b_{1}=1$ then $b_{N}=N-1$ .

e=2+{\cfrac {1}{1+{\cfrac {1}{2+{\cfrac {2}{3+{\cfrac {3}{4+\ddots }}}}}}}}

For Pi, use $a_{0}=3$ then $a_{N}=6$ . $b_{N}=(2N-1)^{2}$ .

\pi =3+{\cfrac {1}{6+{\cfrac {9}{6+{\cfrac {25}{6+\ddots }}}}}}

Ada

<lang Ada>with Ada.Text_IO; use Ada.Text_IO; procedure ContFract is

  type FormType is (Sqrt2, Napier, Pi);
  type Floaty is digits 15;
  package FIO is new Ada.Text_IO.Float_IO (Floaty);

  procedure GetCoefs (form : FormType; n : Natural;
     coefA : out Natural; coefB : out Natural)
  is begin
     case form is
        when Sqrt2 =>
           if n > 0 then coefA := 2; else coefA := 1; end if;
           coefB := 1;
        when Napier =>
           if n > 0 then coefA := n; else coefA := 2; end if;
           if n > 1 then coefB := n - 1; else coefB := 1; end if;
        when Pi =>
           if n > 0 then coefA := 6; else coefA := 3; end if;
           coefB := (2*n - 1)**2;
     end case;
  end GetCoefs;

  function Calc (form : FormType; n : Natural) return Floaty is
     A, B : Natural;
     Temp : Floaty := 0.0;
  begin
     for ni in reverse Natural range 1 .. n loop
        GetCoefs (form, ni, A, B);
        Temp := Floaty (B) / (Floaty (A) + Temp);
     end loop;
     GetCoefs (form, 0, A, B);
     return Floaty (A) + Temp;
  end Calc;

begin

  FIO.Put (Calc (Sqrt2, 200), Exp => 0); New_Line;
  FIO.Put (Calc (Napier, 200), Exp => 0); New_Line;
  FIO.Put (Calc (Pi, 200), Exp => 0); New_Line;

end ContFract; </lang>

Output:

 1.41421356237310
 2.71828182845905
 3.14159262280485

Axiom

Axiom provides a ContinuedFraction domain: <lang Axiom>get(obj) == convergents(obj).1000 -- utility to extract the 1000th value get continuedFraction(1, repeating [1], repeating [2]) :: Float get continuedFraction(2, cons(1,[i for i in 1..]), [i for i in 1..]) :: Float get continuedFraction(3, [i^2 for i in 1.. by 2], repeating [6]) :: Float</lang> Output:<lang Axiom> (1) 1.4142135623 730950488

                                                                 Type: Float

  (2)  2.7182818284 590452354
                                                                 Type: Float

  (3)  3.1415926538 39792926
                                                                 Type: Float</lang>

The value for $\pi$ has an accuracy to only 9 decimal places after 1000 iterations, with an accuracy to 12 decimal places after 10000 iterations.

We could re-implement this, with the same output:<lang Axiom>cf(initial, a, b, n) ==

 n=1 => initial
 temp := 0
 for i in (n-1)..1 by -1 repeat
   temp := a.i/(b.i+temp)
 initial+temp

cf(1, repeating [1], repeating [2], 1000) :: Float cf(2, cons(1,[i for i in 1..]), [i for i in 1..], 1000) :: Float cf(3, [i^2 for i in 1.. by 2], repeating [6], 1000) :: Float</lang>

CoffeeScript

<lang coffeescript># Compute a continuous fraction of the form

a0 + b1 / (a1 + b2 / (a2 + b3 / ...

continuous_fraction = (f) ->

 a = f.a
 b = f.b
 c = 1
 for n in [100000..1]
   c = b(n) / (a(n) + c)
 a(0) + c

A little helper.

p = (a, b) ->

 console.log a
 console.log b
 console.log "---"

do ->

 fsqrt2 =
   a: (n) -> if n is 0 then 1 else 2
   b: (n) -> 1
 p Math.sqrt(2), continuous_fraction(fsqrt2)

 fnapier =
   a: (n) -> if n is 0 then 2 else n
   b: (n) -> if n is 1 then 1 else n - 1
 p Math.E, continuous_fraction(fnapier)

 fpi =
   a: (n) ->
     return 3 if n is 0
     6
   b: (n) ->
     x = 2*n - 1
     x * x
 p Math.PI, continuous_fraction(fpi)

</lang> output

> coffee continued_fraction.coffee 
1.4142135623730951
1.4142135623730951
---
2.718281828459045
2.7182818284590455
---
3.141592653589793
3.141592653589793
---

D

<lang d>import std.typecons;

FP calc(FP, F)(in F fun, in int n) pure nothrow {

   FP temp = 0.0;

   foreach_reverse (ni; 1 .. n+1) {
       immutable a_b = fun(ni);
       temp = cast(FP)a_b[1] / (cast(FP)a_b[0] + temp);
   }
   return cast(FP)fun(0)[0] + temp;

}

// int[2] fsqrt2(in int n) pure nothrow { Tuple!(int,int) fsqrt2(in int n) pure nothrow {

   return tuple(n > 0 ? 2 : 1,
                1);

}

Tuple!(int,int) fnapier(in int n) pure nothrow {

   return tuple(n > 0 ? n : 2,
                n > 1 ? (n - 1) : 1);

}

Tuple!(int,int) fpi(in int n) pure nothrow {

   return tuple(n > 0 ? 6 : 3,
                (2 * n - 1) ^^ 2);

}

void main() {

   import std.stdio;
   writefln("%.19f", calc!real(&fsqrt2, 200));
   writefln("%.19f", calc!real(&fnapier, 200));
   writefln("%.19f", calc!real(&fpi, 200));

}</lang>

Output:

1.4142135623730950487
2.7182818284590452354
3.1415926228048469486

Factor

cfrac-estimate uses rational arithmetic and never truncates the intermediate result. When terms is large, cfrac-estimate runs slow because numerator and denominator grow big.

<lang factor>USING: arrays combinators io kernel locals math math.functions

 math.ranges prettyprint sequences ;

IN: rosetta.cfrac

! Every continued fraction must implement these two words. GENERIC: cfrac-a ( n cfrac -- a ) GENERIC: cfrac-b ( n cfrac -- b )

! square root of 2 SINGLETON: sqrt2 M: sqrt2 cfrac-a

   ! If n is 1, then a_n is 1, else a_n is 2.
   drop { { 1 [ 1 ] } [ drop 2 ] } case ;

M: sqrt2 cfrac-b

   ! Always b_n is 1.
   2drop 1 ;

! Napier's constant SINGLETON: napier M: napier cfrac-a

   ! If n is 1, then a_n is 2, else a_n is n - 1. 
   drop { { 1 [ 2 ] } [ 1 - ] } case ;

M: napier cfrac-b

   ! If n is 1, then b_n is 1, else b_n is n - 1.
   drop { { 1 [ 1 ] } [ 1 - ] } case ;

SINGLETON: pi M: pi cfrac-a

   ! If n is 1, then a_n is 3, else a_n is 6.
   drop { { 1 [ 3 ] } [ drop 6 ] } case ;

M: pi cfrac-b

   ! Always b_n is (n * 2 - 1)^2.
   drop 2 * 1 - 2 ^ ;

cfrac-estimate ( cfrac terms -- number )

   terms cfrac cfrac-a             ! top = last a_n
   terms 1 - 1 [a,b] [ :> n
       n cfrac cfrac-b swap /      ! top = b_n / top
       n cfrac cfrac-a +           ! top = top + a_n
   ] each ;

decimalize ( rational prec -- string )

   rational 1 /mod             ! split whole, fractional parts
   prec 10^ *                  ! multiply fraction by 10 ^ prec
   [ >integer unparse ] bi@    ! convert digits to strings
   :> fraction
   "."                         ! push decimal point
   prec fraction length -
   dup 0 < [ drop 0 ] when
   "0" <repetition> concat     ! push padding zeros
   fraction 4array concat ;

<PRIVATE

main ( -- )

   " Square root of 2: " write
   sqrt2 50 cfrac-estimate 30 decimalize print
   "Napier's constant: " write
   napier 50 cfrac-estimate 30 decimalize print
   "               Pi: " write
   pi 950 cfrac-estimate 10 decimalize print ;

PRIVATE>

MAIN: main</lang>

Output:

 Square root of 2: 1.414213562373095048801688724209
Napier's constant: 2.718281828459045235360287471352
               Pi: 3.1415926538

Forth

Translation of: D

<lang forth>: fsqrt2 1 s>f 0> if 2 s>f else fdup then ;

fnapier dup dup 1 > if 1- else drop 1 then s>f dup 1 < if drop 2 then s>f ;

fpi dup 2* 1- dup * s>f 0> if 6 else 3 then s>f ;

                                      ( n -- f1 f2)

cont.fraction ( xt n -- f)

 1 swap 1+ 0 s>f                      \ calculate for 1 .. n
 do i over execute frot f+ f/ -1 +loop
 0 swap execute fnip f+               \ calcucate for 0

</lang>

Output:

' fsqrt2  200 cont.fraction f. cr 1.4142135623731
 ok
' fnapier 200 cont.fraction f. cr 2.71828182845905
 ok
' fpi     200 cont.fraction f. cr 3.14159268391981
 ok

Go

<lang go>package main

import "fmt"

type cfTerm struct {

   a, b int

}

// follows subscript convention of mathworld and WP where there is no b(0). // cf[0].b is unused in this representation. type cf []cfTerm

func cfSqrt2(nTerms int) cf {

   f := make(cf, nTerms)
   for n := range f {
       f[n] = cfTerm{2, 1}
   }
   f[0].a = 1
   return f

}

func cfNap(nTerms int) cf {

   f := make(cf, nTerms)
   for n := range f {
       f[n] = cfTerm{n, n - 1}
   }
   f[0].a = 2
   f[1].b = 1
   return f

}

func cfPi(nTerms int) cf {

   f := make(cf, nTerms)
   for n := range f {
       g := 2*n - 1
       f[n] = cfTerm{6, g * g}
   }
   f[0].a = 3
   return f

}

func (f cf) real() (r float64) {

   for n := len(f) - 1; n > 0; n-- {
       r = float64(f[n].b) / (float64(f[n].a) + r)
   }
   return r + float64(f[0].a)

}

func main() {

   fmt.Println("sqrt2:", cfSqrt2(20).real())
   fmt.Println("nap:  ", cfNap(20).real())
   fmt.Println("pi:   ", cfPi(20).real())

}</lang>

Output:

sqrt2: 1.4142135623730965
nap:   2.7182818284590455
pi:    3.141623806667839

Haskell

<lang haskell>import Data.List (unfoldr) import Data.Char (intToDigit)

-- continued fraction represented as a (possibly infinite) list of pairs sqrt2, napier, myPi :: [(Integer, Integer)] sqrt2 = zip (1 : [2,2..]) [1,1..] napier = zip (2 : [1..]) (1 : [1..]) myPi = zip (3 : [6,6..]) (map (^2) [1,3..])

-- approximate a continued fraction after certain number of iterations approxCF :: (Integral a, Fractional b) => Int -> [(a, a)] -> b approxCF t =

 foldr (\(a,b) z -> fromIntegral a + fromIntegral b / z) 1 . take t

-- infinite decimal representation of a real number decString :: RealFrac a => a -> String decString frac = show i ++ '.' : decString' f where

 (i,f) = properFraction frac
 decString' = map intToDigit . unfoldr (Just . properFraction . (10*))

main :: IO () main = mapM_ (putStrLn . take 200 . decString .

             (approxCF 950 :: [(Integer, Integer)] -> Rational))
            [sqrt2, napier, myPi]</lang>

Output:

1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572735013846230912297024924836055850737212644121497099935831413222665927505592755799950501152782060571
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019
3.141592653297590947683406834261190738869139611505752231394089152890909495973464508817163306557131591579057202097715021166512662872910519439747609829479577279606075707015622200744006783543589980682386

J

  cfrac=: +`% / NB. Evaluate a list as a continued fraction

  sqrt2=: cfrac 1 1,200$2 1x
  pi=:cfrac 3, , ,&6"0 *:<:+:>:i.100x
  e=: cfrac 2 1, , ,~"0 >:i.100x

  NB. translate from fraction to decimal string
  NB. translated from factor
  dec =: (-@:[ (}.,'.',{.) ":@:<.@:(* 10x&^)~)"0

  100 10 100 dec sqrt2, pi, e

1.4142135623730950488016887242096980785696718753769480731766797379907324784621205551109457595775322165 3.1415924109 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274 </lang>

Maxima

<lang maxima>cfeval(x) := block([a, b, n, z], a: x[1], b: x[2], n: length(a), z: 0,

  for i from n step -1 thru 2 do z: b[i]/(a[i] + z), a[1] + z)$

cf_sqrt2(n) := [cons(1, makelist(2, i, 2, n)), cons(0, makelist(1, i, 2, n))]$

cf_e(n) := [cons(2, makelist(i, i, 1, n - 1)), append([0, 1], makelist(i, i, 1, n - 2))]$

cf_pi(n) := [cons(3, makelist(6, i, 2, n)), cons(0, makelist((2*i - 1)^2, i, 1, n - 1))]$

cfeval(cf_sqrt2(20)), numer; /* 1.414213562373097 */ % - sqrt(2), numer; /* 1.3322676295501878*10^-15 */

cfeval(cf_e(20)), numer; /* 2.718281828459046 */ % - %e, numer; /* 4.4408920985006262*10^-16 */

cfeval(cf_pi(20)), numer; /* 3.141623806667839 */ % - %pi, numer; /* 3.115307804568701*10^-5 */

/* convergence is much slower for pi */ fpprec: 20$ x: cfeval(cf_pi(10000))$ bfloat(x - %pi); /* 2.4999999900104930006b-13 */</lang>

PARI/GP

Partial solution for simple continued fractions. <lang parigp>back(v)=my(t=contfracpnqn(v));t[1,1]/t[2,1]*1. back(vector(100,i,2-(i==1)))</lang>

Output:

%1 = 1.4142135623730950488016887242096980786

PL/I

<lang PL/I>/* Version for SQRT(2) */ test: proc options (main);

  declare n fixed;

denom: procedure (n) recursive returns (float (18));

  declare n fixed;
  n = n + 1;
  if n > 100 then return (2);
  return (2 + 1/denom(n));

end denom;

  put (1 + 1/denom(2));

end test; </lang> Result:

 1.41421356237309505E+0000

Version for SQRT2, NAPIER, PI <lang>/* Derived from continued fraction in Wiki Ada program */

continued_fractions: /* 6 Sept. 2012 */

  procedure options (main);
  declare (Sqrt2 initial (1), napier initial (2), pi initial (3)) fixed (1);

Get_Coeffs: procedure (form, n, coefA, coefB);

     declare form fixed (1), n fixed, (coefA, coefB) float (18);

     select (form);
        when (Sqrt2) do;
              if n > 0 then coefA = 2; else coefA = 1;
              coefB = 1;
           end;
        when (Napier) do;
              if n > 0 then coefA = n; else coefA = 2;
              if n > 1 then coefB = n - 1; else coefB = 1;
           end;
        when (Pi) do;
              if n > 0 then coefA = 6; else coefA = 3;
              coefB = (2*n - 1)**2;
           end;
     end;
  end Get_Coeffs;

  Calc: procedure (form, n) returns (float (18));
     declare form fixed (1), n fixed;
     declare (A, B) float (18);
     declare Temp float (18) initial (0);
     declare ni fixed;

     do ni = n to 1 by -1;
        call Get_Coeffs (form, ni, A, B);
        Temp = B/(A + Temp);
     end;
     call Get_Coeffs (form, 0, A, B);
     return (A + Temp);
  end Calc;

  put      edit ('SQRT2=',  calc(sqrt2,  200)) (a(10), f(20,17));
  put skip edit ('NAPIER=', calc(napier, 200)) (a(10), f(20,17));
  put skip edit ('PI=',     calc(pi,     200)) (a(10), f(20,17));

end continued_fractions; </lang> Output:

SQRT2=     1.41421356237309505
NAPIER=    2.71828182845904524
PI=        3.14159262280484695

Prolog

<lang Prolog>continued_fraction :- % square root 2 continued_fraction(200, sqrt_2_ab, V1), format('sqrt(2) = ~w~n', [V1]),

% napier continued_fraction(200, napier_ab, V2), format('e = ~w~n', [V2]),

% pi continued_fraction(200, pi_ab, V3), format('pi = ~w~n', [V3]).

% code for continued fractions continued_fraction(N, Compute_ab, V) :- continued_fraction(N, Compute_ab, 0, V).

continued_fraction(0, Compute_ab, Temp, V) :- call(Compute_ab, 0, A, _), V is A + Temp.

continued_fraction(N, Compute_ab, Tmp, V) :- call(Compute_ab, N, A, B), Tmp1 is B / (A + Tmp), N1 is N - 1, continued_fraction(N1, Compute_ab, Tmp1, V).

% specific codes for examples % definitions for square root of 2 sqrt_2_ab(0, 1, 1). sqrt_2_ab(_, 2, 1).

% definitions for napier napier_ab(0, 2, _). napier_ab(1, 1, 1). napier_ab(N, N, V) :- V is N - 1.

% definitions for pi pi_ab(0, 3, _). pi_ab(N, 6, V) :- V is (2 * N - 1)*(2 * N - 1). </lang> Output :

 ?- continued_fraction.
sqrt(2) = 1.4142135623730951
e       = 2.7182818284590455
pi      = 3.141592622804847
true .

Python

Works with: Python version 2.6+ and 3.x

<lang python> from fractions import Fraction import itertools try: zip = itertools.izip except: pass

The Continued Fraction

def CF(a, b, t):

 terms = list(itertools.islice(zip(a, b), t))
 z = Fraction(1,1)
 for a, b in reversed(terms):
   z = a + b / z
 return z

Approximates a fraction to a string

def pRes(x, d):

 q, x = divmod(x, 1)
 res = str(q)
 res += "."
 for i in range(d):
   x *= 10
   q, x = divmod(x, 1)
   res += str(q)
 return res

Test the Continued Fraction for sqrt2

def sqrt2_a():

 yield 1
 for x in itertools.repeat(2):
   yield x

def sqrt2_b():

 for x in itertools.repeat(1):
   yield x

cf = CF(sqrt2_a(), sqrt2_b(), 950) print(pRes(cf, 200))

1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273501384623091229702492483605585073721264412149709993583141322266592750559275579995050115278206057147

Test the Continued Fraction for Napier's Constant

def Napier_a():

 yield 2
 for x in itertools.count(1):
   yield x

def Napier_b():

 yield 1
 for x in itertools.count(1):
   yield x

cf = CF(Napier_a(), Napier_b(), 950) print(pRes(cf, 200))

2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526059563073813232862794349076323382988075319525101901

Test the Continued Fraction for Pi

def Pi_a():

 yield 3
 for x in itertools.repeat(6):
   yield x

def Pi_b():

 for x in itertools.count(1,2):
   yield x*x

cf = CF(Pi_a(), Pi_b(), 950) print(pRes(cf, 10))

3.1415926532

</lang>

Fast iterative version

Translation of: D

<lang python>from decimal import Decimal, getcontext

def calc(fun, n):

   temp = Decimal("0.0")

   for ni in xrange(n+1, 0, -1):
       (a, b) = fun(ni)
       temp = Decimal(b) / (a + temp)

   return fun(0)[0] + temp

def fsqrt2(n):

   return (2 if n > 0 else 1, 1)

def fnapier(n):

   return (n if n > 0 else 2, (n - 1) if n > 1 else 1)

def fpi(n):

   return (6 if n > 0 else 3, (2 * n - 1) ** 2)

getcontext().prec = 50 print calc(fsqrt2, 200) print calc(fnapier, 200) print calc(fpi, 200)</lang>

Output:

1.4142135623730950488016887242096980785696718753770
2.7182818284590452353602874713526624977572470937000
3.1415926839198062649342019294083175420335002640134

REXX

The $CF subroutine (for continued fractions) isn't limited to positive integers,
any form of REXX numbers (negative, exponentiated, decimal fractions) can be used.

There isn't any practical limit on the precision that can be used, although 100k digits would be a bit unwieldly to display. <lang rexx>/*REXX program calculates & displays values of some continued fractions.*/ numeric digits 100 /*use a hundred digits for displays. */ terms=500 /*use 5 hundred terms for calculations*/ /*══════════════════════════════════════════════════════════════════════*/ a=copies(' 2',terms); call tell '√2',$cf(1 a) /*══════════════════════════════════════════════════════════════════════*/

   do j=1 for terms;  a=a j;           end;  call tell 'e',$cf(2 a,1 a)

/*══════════════════════════════════════════════════════════════════════*/ a=copies(' 1',terms); call tell 'φ, phi',$cf(1 a) /*══════════════════════════════════════════════════════════════════════*/ a=copies(' 1 2',terms%2); call tell '√3',$cf(1 a) /*also: 2∙sin(π/3) */ /*══════════════════════════════════════════════════════════════════════*/

   do j=1 for terms%2 by 2; a=a j 1;   end;  call tell 'tan(1)',$cf(1 a)

/*══════════════════════════════════════════════════════════════════════*/

   do j=1 for terms;        a=a 2*j+1; end;  call tell 'coth(1)',$cf(1 a)

/*══════════════════════════════════════════════════════════════════════*/

   do j=1 for terms;        a=a 4*j+2; end;  call tell 'coth(½)',$cf(2 a)     /*also:  [e+1] ÷ [e-1] */

/*══════════════════════════════════════════════════════════════════════*/ terms=100000 /*use 100,000 terms for π calc.*/ a=copies(' 6',terms)

   do j=1 for terms;  b=b (2*j-1)**2;  end;  call tell 'π, pi',$cf(3 a,b)

exit /*stick a fork in it, we're done.*/ /*────────────────────────────────$CF subroutine────────────────────────*/ $cf: procedure; parse arg C x,y; !=0; numeric digits digits()*2

          do k=words(x) to 1 by -1;   a=word(x,k);  b=word(word(y,k) 1,1)
          d=a+!;  if d=0 then call divZero  /*in case divisor is bogus.*/
          !=b/d
          end   /*k*/

return !+C /*────────────────────────────────DIVZERO subroutine────────────────────*/ divZero: say; say '***error!***'; say 'division by zero.'; say; exit 13 /*────────────────────────────────TELL subroutine────────────────────────────────*/ tell:parse arg _,v;say right(_,7) '=' format(v) ' α terms='left(a,50);a=;b=;return</lang> output

     √2 = 1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573   α terms= 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
      e = 2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427   α terms= 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 2
 φ, phi = 1.618033988749894848204586834365638117720309179805762862135448622705260462818902449707207204189391137   α terms= 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
     √3 = 1.732050807568877293527446341505872366942805253810380628055806979451933016908800037081146186757248576   α terms= 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1
 tan(1) = 1.557407724654902230506974807458360173087250772381520038383946605698861397151727289555099965202242984   α terms= 1 1 3 1 5 1 7 1 9 1 11 1 13 1 15 1 17 1 19 1 21 1
coth(1) = 1.313035285499331303636161246930847832912013941240452655543152967567084270461874382674679241480856303   α terms= 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37
coth(½) = 2.163953413738652848770004010218023117093738602150792272533574119296087634783339486574409418809750115   α terms= 6 10 14 18 22 26 30 34 38 42 46 50 54 58 62 66 70
  π, pi = 3.141592653589792988470143264530440384041017830472772036746332303472711537960073664096818977224037083   α terms= 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6

Note: even with 400 digit accuracy and 100,000 terms, the calculate of π (pi) is only accurate to 15 digits.

Ruby

<lang ruby>require 'bigdecimal'

square root of 2

sqrt2 = Object.new def sqrt2.a(n); n == 1 ? 1 : 2; end def sqrt2.b(n); 1; end

Napier's constant

napier = Object.new def napier.a(n); n == 1 ? 2 : n - 1; end def napier.b(n); n == 1 ? 1 : n - 1; end

pi = Object.new def pi.a(n); n == 1 ? 3 : 6; end def pi.b(n); (2*n - 1)**2; end

Estimates the value of a continued fraction _cfrac_, to _prec_
decimal digits of precision. Returns a BigDecimal. _cfrac_ must
respond to _cfrac.a(n)_ and _cfrac.b(n)_ for integer _n_ >= 1.

def estimate(cfrac, prec)

 last_result = nil
 terms = prec

 loop do
   # Estimate continued fraction for _n_ from 1 to _terms_.
   result = cfrac.a(terms)
   (terms - 1).downto(1) do |n|
     a = BigDecimal cfrac.a(n)
     b = BigDecimal cfrac.b(n)
     digits = [b.div(result, 1).exponent + prec, 1].max
     result = a + b.div(result, digits)
   end
   result = result.round(prec)

   if result == last_result
     return result
   else
     # Double _terms_ and try again.
     last_result = result
     terms *= 2
   end
 end

end

puts estimate(sqrt2, 50).to_s('F') puts estimate(napier, 50).to_s('F') puts estimate(pi, 10).to_s('F')</lang>

Output:

$ ruby cfrac.rb                                                              
1.41421356237309504880168872420969807856967187537695
2.71828182845904523536028747135266249775724709369996
3.1415926536

Scala

Works with: Scala version 2.9.1

Note that Scala-BigDecimal provides a precision of 34 digits. Therefore we take a limitation of 32 digits to avoiding rounding problems. <lang Scala>object CF extends App {

 import Stream._
 val sqrt2 = 1 #:: from(2,0) zip from(1,0)
 val napier = 2 #:: from(1) zip (1 #:: from(1))
 val pi = 3 #:: from(6,0) zip (from(1,2) map {x=>x*x})

 // reference values, source: wikipedia
 val refPi     = "3.14159265358979323846264338327950288419716939937510"
 val refNapier = "2.71828182845904523536028747135266249775724709369995"
 val refSQRT2  = "1.41421356237309504880168872420969807856967187537694"

 def calc(cf: Stream[(Int, Int)], numberOfIters: Int=200): BigDecimal = {
   (cf take numberOfIters toList).foldRight[BigDecimal](1)((a, z) => a._1+a._2/z)
 }
 
 def approx(cfV: BigDecimal, cfRefV: String): String = {
   val p: Pair[Char,Char] => Boolean = pair =>(pair._1==pair._2)
   ((cfV.toString+" "*34).substring(0,34) zip cfRefV.toString.substring(0,34))
     .takeWhile(p).foldRight[String]("")((a:Pair[Char,Char],z)=>a._1+z)
 }

 List(("sqrt2",sqrt2,50,refSQRT2),("napier",napier,50,refNapier),("pi",pi,3000,refPi)) foreach {t=>
   val (name,cf,iters,refV) = t
   val cfV = calc(cf,iters)
   println(name+":")
   println("ref value: "+refV.substring(0,34))
   println("cf value:  "+(cfV.toString+" "*34).substring(0,34))
   println("precision: "+approx(cfV,refV))
   println()
 }

}</lang> Output:

sqrt2:
ref value: 1.41421356237309504880168872420969
cf value:  1.41421356237309504880168872420969
precision: 1.41421356237309504880168872420969

napier:
ref value: 2.71828182845904523536028747135266
cf value:  2.71828182845904523536028747135266
precision: 2.71828182845904523536028747135266

pi:
ref value: 3.14159265358979323846264338327950
cf value:  3.14159265358052780404906362935452
precision: 3.14159265358

For higher accuracy of pi we have to take more iterations. Unfortunately the foldRight function in calc isn't tail recursiv - therefore a stack overflow exception will be thrown for higher numbers of iteration, thus we have to implement an iterative way for calculation: <lang Scala>object CFI extends App {

 import Stream._
 val sqrt2 = 1 #:: from(2,0) zip from(1,0)
 val napier = 2 #:: from(1) zip (1 #:: from(1))
 val pi = 3 #:: from(6,0) zip (from(1,2) map {x=>x*x})

 // reference values, source: wikipedia
 val refPi     = "3.14159265358979323846264338327950288419716939937510"
 val refNapier = "2.71828182845904523536028747135266249775724709369995"
 val refSQRT2  = "1.41421356237309504880168872420969807856967187537694"

 def calc_i(cf: Stream[(Int, Int)], numberOfIters: Int=50): BigDecimal = {
   val cfl = cf take numberOfIters toList
   var z: BigDecimal = 1.0
   for (i <- 0 to cfl.size-1 reverse) 
     z=cfl(i)._1+cfl(i)._2/z
   z
 }
 
 def approx(cfV: BigDecimal, cfRefV: String): String = {
   val p: Pair[Char,Char] => Boolean = pair =>(pair._1==pair._2)
   ((cfV.toString+" "*34).substring(0,34) zip cfRefV.toString.substring(0,34))
     .takeWhile(p).foldRight[String]("")((a:Pair[Char,Char],z)=>a._1+z)
 }

 List(("sqrt2",sqrt2,50,refSQRT2),("napier",napier,50,refNapier),("pi",pi,50000,refPi)) foreach {t=>
   val (name,cf,iters,refV) = t
   val cfV = calc_i(cf,iters)
   println(name+":")
   println("ref value: "+refV.substring(0,34))
   println("cf value:  "+(cfV.toString+" "*34).substring(0,34))
   println("precision: "+approx(cfV,refV))
   println()
 }

}</lang> Output:

sqrt2:
ref value: 1.41421356237309504880168872420969
cf value:  1.41421356237309504880168872420969
precision: 1.41421356237309504880168872420969

napier:
ref value: 2.71828182845904523536028747135266
cf value:  2.71828182845904523536028747135266
precision: 2.71828182845904523536028747135266

pi:
ref value: 3.14159265358979323846264338327950
cf value:  3.14159265358983426214354599901745
precision: 3.141592653589

Tcl

Works with: tclsh version 8.6

Translation of: Python

Note that Tcl does not provide arbitrary precision floating point numbers by default, so all result computations are done with IEEE doubles. <lang tcl>package require Tcl 8.6

Term generators; yield list of pairs

proc r2 {} {

   yield {1 1}
   while 1 {yield {2 1}}

} proc e {} {

   yield {2 1}
   while 1 {yield [list [incr n] $n]}

} proc pi {} {

   set n 0; set a 3
   while 1 {

yield [list $a [expr {(2*[incr n]-1)**2}]] set a 6

}

Continued fraction calculator

proc cf {generator {termCount 50}} {

   # Get the chunk of terms we want to work with
   set terms [list [coroutine cf.c $generator]]
   while {[llength $terms] < $termCount} {

lappend terms [cf.c]

   }
   rename cf.c {}

   # Merge the terms to compute the result
   set val 0.0
   foreach pair [lreverse $terms] {

lassign $pair a b set val [expr {$a + $b/$val}]

   }
   return $val

}

Demonstration

puts [cf r2] puts [cf e] puts [cf pi 250]; # Converges more slowly</lang>

Output:

1.4142135623730951
2.7182818284590455
3.1415926373965735

XPL0

The number of iterations (N) needed to get the 13 digits of accuracy was determined by experiment. <lang XPL0>include c:\cxpl\codes; int N; real A, B, F; [Format(1, 15); A:= 2.0; B:= 1.0; N:= 16; IntOut(0, N); CrLf(0); F:= 0.0; while N>=1 do [F:= B/(A+F); N:= N-1]; RlOut(0, 1.0+F); CrLf(0); RlOut(0, sqrt(2.0)); CrLf(0);

N:= 13; IntOut(0, N); CrLf(0); F:= 0.0; while N>=2 do [F:= float(N-1)/(float(N)+F); N:= N-1]; RlOut(0, 2.0 + 1.0/(1.0+F)); CrLf(0); RlOut(0, Exp(1.0)); CrLf(0);

N:= 10000; IntOut(0, N); CrLf(0); F:= 0.0; while N>=1 do [F:= float(sq(2*N-1))/(6.0+F); N:= N-1]; RlOut(0, 3.0+F); CrLf(0); RlOut(0, ACos(-1.0)); CrLf(0); ]</lang>

Output:

16
1.414213562372820
1.414213562373100
13
2.718281828459380
2.718281828459050
10000
3.141592653589540
3.141592653589790