Continued fraction

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, 10000), Exp => 0); New_Line;

end ContFract; </lang>

Output:

 1.41421356237310
 2.71828182845905

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>

C++

<lang cpp>#include <iomanip>

include <iostream>
include <tuple>

typedef std::tuple<double,double> coeff_t; // coefficients type typedef coeff_t (*func_t)(int); // callback function type

double calc(func_t func, int n) {

   double a, b, temp = 0;
   for (; n > 0; --n) {
       std::tie(a, b) = func(n);
       temp = b / (a + temp);
   }
   std::tie(a, b) = func(0);
   return a + temp;

}

coeff_t sqrt2(int n) {

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

}

coeff_t napier(int n) {

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

}

coeff_t pi(int n) {

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

}

int main() {

   std::streamsize old_prec = std::cout.precision(15); // set output digits
   std::cout 
       << calc(sqrt2, 20) << '\n'
       << calc(napier, 15) << '\n'
       << calc(pi, 10000) << '\n'
       << std::setprecision(old_prec); // reset precision

}</lang>

Output:

1.41421356237309
2.71828182845905
3.14159265358954

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;

alias Tuple!(int,"a", int,"b") Pair;

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

   FP temp = 0;

   foreach_reverse (ni; 1 .. n+1) {
       immutable p = fun(ni);
       temp = cast(FP)p.b / (cast(FP)p.a + temp);
   }
   return cast(FP)fun(0).a + temp;

}

// int[2] fsqrt2(in int n) pure nothrow { Pair fsqrt2(in int n) pure nothrow {

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

}

Pair fnapier(in int n) pure nothrow {

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

}

Pair fpi(in int n) pure nothrow {

   return Pair(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

<lang haskell>import Data.Ratio

-- ignoring the task-given pi sequence: sucky convergence -- pie = zip (3:repeat 6) (map (^2) [1,3..])

pie = zip (0:[1,3..]) (4:map (^2) [1..]) sqrt2 = zip (1:repeat 2) (repeat 1) napier = zip (2:[1..]) (1:[1..])

-- truncate after n terms cf2rat n = foldr (\(a,b) f -> (a%1) + ((b%1) / f)) (1%1) . take n

-- truncate after error is at most 1/p cf2rat_p p s = f $ map (\i -> (cf2rat i s, cf2rat (1+i) s)) $ map (2^) [0..] where f ((x,y):ys) = if abs (x-y) < (1/fromIntegral p) then x else f ys

-- returns a decimal string of n digits after the dot; all digits should -- be correct (doesn't mean it's the best approximation! the decimal -- string is simply truncated to given digits: pi=3.141 instead of 3.142) cf2dec n = (ratstr n) . cf2rat_p (10^n) where ratstr l a = (show t) ++ '.':fracstr l n d where d = denominator a (t, n) = quotRem (numerator a) d fracstr 0 _ _ = [] fracstr l n d = (show t)++ fracstr (l-1) n1 d where (t,n1) = quotRem (10 * n) d

main = do putStrLn $ cf2dec 200 sqrt2 putStrLn $ cf2dec 200 napier putStrLn $ cf2dec 200 pie</lang>

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>

NetRexx

<lang netrexx>/* REXX ***************************************************************

Derived from REXX ... Derived from PL/I with a little "massage"
SQRT2= 1.41421356237309505 <- PL/I Result
1.41421356237309504880168872421 <- NetRexx Result 30 digits
NAPIER= 2.71828182845904524
2.71828182845904523536028747135
PI= 3.14159262280484695
3.14159262280484694855146925223
07.09.2012 Walter Pachl
08.09.2012 Walter Pachl simplified (with the help of a friend)
- - - - /

options replace format comments java crossref savelog symbols

 class CFB public

properties static

 Numeric Digits 30
 Sqrt2 =1
 napier=2
 pi    =3
 a     =0
 b     =0

method main(args = String[]) public static

 Say 'SQRT2='.left(7)  calc(sqrt2,  200)
 Say 'NAPIER='.left(7) calc(napier, 200)
 Say 'PI='.left(7)     calc(pi,     200)
 Return

method get_Coeffs(form,n) public static

 select
   when form=Sqrt2 Then do
     if n > 0 then a = 2; else a = 1
     b = 1
     end
   when form=Napier Then do
     if n > 0 then a = n; else a = 2
     if n > 1 then b = n - 1; else b = 1
     end
   when form=pi Then do
     if n > 0 then a = 6; else a = 3
     b = (2*n - 1)**2
     end
   end
 Return

method calc(form,n) public static

 temp=0
 loop ni = n to 1 by -1
   Get_Coeffs(form,ni)
   temp = b/(a + temp)
   end
 Get_Coeffs(form,0)
 return (a + temp)</lang>

Who could help me make a,b,sqrt2,napier,pi global (public) variables? This would simplify the solution:-)

I got this help and simplified the program.

However, I am told that 'my' value of pi is incorrect. I will investigate!

Apparently the coefficients given in the task description are only good for an approxunation. One should, therefore, not SHOW more that 15 digits. See http://de.wikipedia.org/wiki/Kreiszahl

See Rexx for a better computation

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 NAPIER: <lang>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 (n);
  return (n + n/denom(n));

end denom;

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

end test;</lang>

 2.71828182845904524E+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,   99999)) (a(10), f(20,17));

end continued_fractions; </lang> Output:

SQRT2=     1.41421356237309505
NAPIER=    2.71828182845904524
PI=        3.14159265358979349

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

Version 1

The CF subroutine (for continued fractions) isn't limited to positive integers.
Any form of REXX numbers (negative, exponentiated, decimal fractions) can be used;
[note the use of negative fractions for the ß terms when computing √½].

There isn't any practical limit on the precision that can be used, although 100k digits would be a bit unwieldly to display.

A generalized √ function was added to calculate a few low integers (and also ½).
In addition, ½π was calculated (as described in the talk page under Gold Credit).

More code is used for nicely formatting the output than the continued fraction calculation. <lang rexx>/*REXX program calculates and displays values of some specific continued*/ /*───────────── fractions (along with their α and ß terms). */ /*───────────── Continued fractions: also known as anthyphairetic ratio.*/ T=500 /*use 500 terms for calculations.*/ showDig=100; numeric digits 2*showDig /*use 100 digits for the display.*/ a=; @=; b= /*omitted ß terms are assumed = 1*/ /*══════════════════════════════════════════════════════════════════════*/ a=1 rep(2); call tell '√2' /*══════════════════════════════════════════════════════════════════════*/ a=1 rep(1 2); call tell '√3' /*also: 2∙sin(π/3) */ /*══════════════════════════════════════ ___ ════════════════════════*/

                          /*generalized  √ N   */
     do N=2 to 11;        a=1 rep(2); b=rep(N-1);      call tell 'gen √'N;   end

N=1/2; a=1 rep(2); b=rep(N-1); call tell 'gen √½' /*══════════════════════════════════════════════════════════════════════*/

     do j=1 for T;        a=a j; end;   b=1 a; a=2 a;  call tell 'e'

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

     do j=1 for T by 2;   a=a j;        b=b j+1; end;  call tell '1÷[√e-1]'

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

     do j=1 for T;        a=a j;     end; b=a; a=0 a;  call tell '1÷[e-1]'

/*══════════════════════════════════════════════════════════════════════*/ a=1 rep(1); call tell 'φ, phi' /*══════════════════════════════════════════════════════════════════════*/ a=1; do j=1 for T by 2; a=a j 1; end; call tell 'tan(1)' /*══════════════════════════════════════════════════════════════════════*/ a=1; do j=1 for T; a=a 2*j+1; end; call tell 'coth(1)' /*══════════════════════════════════════════════════════════════════════*/ a=2; do j=1 for T; a=a 4*j+2; end; call tell 'coth(½)' /*also: [e+1] ÷ [e-1] */ /*══════════════════════════════════════════════════════════════════════*/ T=10000 a=1 rep(2)

     do j=1 for T by 2;   b=b j**2;  end;              call tell '4÷π'

/*══════════════════════════════════════════════════════════════════════*/ T=10000 a=1; do j=1 for T; a=a 1/j; @=@ '1/'j; end; call tell '½π, ½pi' /*══════════════════════════════════════════════════════════════════════*/ T=10000 a=0 1 rep(2)

     do j=1 for T by 2;   b=b j**2;     end;   b=4 b;  call tell 'π, pi'

/*══════════════════════════════════════════════════════════════════════*/ T=10000 a=0; do j=1 for T; a=a j*2-1; b=b j**2; end; b=4 b; call tell 'π, pi' /*══════════════════════════════════════════════════════════════════════*/ T=100000 a=3 rep(6)

     do j=1 for T by 2;       b=b j**2; end;           call tell 'π, pi'

exit /*stick a fork in it, we're done.*/

/*────────────────────────────────CF subroutine─────────────────────────*/ cf: procedure; parse arg C x,y; !=0; numeric digits digits()+5

          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                             /*here's a binary mosh pit.*/
          end   /*k*/

return !+C /*────────────────────────────────DIVZERO subroutine────────────────────*/ divZero: say; say '***error!***'; say 'division by zero.'; say; exit 13 /*────────────────────────────────GETT subroutine───────────────────────*/ getT: parse arg stuff,width,ma,mb,_

               do m=1; mm=m+ma; mn=max(1,m-mb);   w=word(stuff,m)
               w=right(w,max(length(word(As,mm)),length(word(Bs,mn)),length(w)))
               if length(_ w)>width  then leave   /*stop getting terms?*/
               _=_ w                              /*whole, don't chop. */
               end   /*m*/                        /*done building terms*/

return strip(_) /*strip leading blank*/ /*────────────────────────────────REP subroutine────────────────────────*/ rep: parse arg rep; return space(copies(' 'rep, T%words(rep))) /*────────────────────────────────RF subroutine─────────────────────────*/ rf: parse arg xxx,z

        do m=1 for T; w=word(xxx,m)    ;  if w=='1/1' | w=1  then w=1
        if w=='1/2' | w=1/2  then w='½';  if w=-.5           then w='-½'
        if w=='1/4' | w=1/4  then w='¼';  if w=-.24          then w='-¼'
        z=z w
        end

return z /*done re-formatting.*/ /*────────────────────────────────TELL subroutine───────────────────────*/ tell: parse arg ?; v=cf(a,b); numeric digits showdig; As=rf(@ a); Bs=rf(b)

     say right(?,8)   '='   left(v/1,showdig)   '  α terms= '   getT(As,72  ,0,1)
     if b\== then say right(,8+2+showdig+1) '  ß terms=   ' getT(Bs,72-2,1,0)
     a=;  @=;  b=;   return</lang>

output

      √2 = 1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157   α terms=  1 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 2 2 2 2 2 2 2 2 2 2
      √3 = 1.73205080756887729352744634150587236694280525381038062805580697945193301690880003708114618675724857   α terms=  1 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 2 1 2 1 2 1 2 1 2 1
  gen √2 = 1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157   α terms=  1 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 2 2 2 2 2 2 2 2 2 2
                                                                                                                  ß 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 1 1 1 1 1 1 1 1 1 1
  gen √3 = 1.73205080756887729352744634150587236694280525381038062805580697945193301690880003708114618675724857   α terms=  1 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 2 2 2 2 2 2 2 2 2 2
                                                                                                                  ß 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 2 2 2 2 2 2 2 2 2 2
  gen √4 = 2                                                                                                      α terms=  1 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 2 2 2 2 2 2 2 2 2 2
                                                                                                                  ß terms=    3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
  gen √5 = 2.23606797749978969640917366873127623544061835961152572427089724541052092563780489941441440837878227   α terms=  1 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 2 2 2 2 2 2 2 2 2 2
                                                                                                                  ß terms=    4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
  gen √6 = 2.44948974278317809819728407470589139196594748065667012843269256725096037745731502653985943310464023   α terms=  1 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 2 2 2 2 2 2 2 2 2 2
                                                                                                                  ß terms=    5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
  gen √7 = 2.64575131106459059050161575363926042571025918308245018036833445920106882323028362776039288647454361   α terms=  1 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 2 2 2 2 2 2 2 2 2 2
                                                                                                                  ß 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 6 6 6 6 6 6 6 6 6 6
  gen √8 = 2.82842712474619009760337744841939615713934375075389614635335947598146495692421407770077506865528314   α terms=  1 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 2 2 2 2 2 2 2 2 2 2
                                                                                                                  ß terms=    7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
  gen √9 = 3                                                                                                      α terms=  1 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 2 2 2 2 2 2 2 2 2 2
                                                                                                                  ß terms=    8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8
 gen √10 = 3.16227766016837933199889354443271853371955513932521682685750485279259443863923822134424810837930029   α terms=  1 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 2 2 2 2 2 2 2 2 2 2
                                                                                                                  ß terms=    9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
 gen √11 = 3.31662479035539984911493273667068668392708854558935359705868214611648464260904384670884339912829065   α terms=  1  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
                                                                                                                  ß terms=    10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
  gen √½ = 0.70710678118654752440084436210484903928483593768847403658833986899536623923105351942519376716382078   α terms=  1  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
                                                                                                                  ß terms=    -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½
       e = 2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642   α terms=  2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
                                                                                                                  ß terms=    1 1 2 3 4 5 6 7 8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
1÷[√e-1] = 1.54149408253679828413110344447251463834045923684188210947413695663754263914331480707182572408500774   α terms=  1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49
                                                                                                                  ß terms=    2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
 1÷[e-1] = 0.58197670686932642438500200510901155854686930107539613626678705964804381739166974328720470940487505   α terms=  0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
                                                                                                                  ß terms=    1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
  φ, phi = 1.61803398874989484820458683436563811772030917980576286213544862270526046281890244970720720418939113   α 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 1 1 1 1 1 1 1 1 1 1 1
  tan(1) = 1.55740772465490223050697480745836017308725077238152003838394660569886139715172728955509996520224298   α terms=  1 1 1 3 1 5 1 7 1 9 1 11 1 13 1 15 1 17 1 19 1 21 1 23 1 25 1 27 1 29 1
 coth(1) = 1.31303528549933130363616124693084783291201394124045265554315296756708427046187438267467924148085630   α terms=  1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49
 coth(½) = 2.16395341373865284877000401021802311709373860215079227253357411929608763478333948657440941880975011   α terms=  2 6 10 14 18 22 26 30 34 38 42 46 50 54 58 62 66 70 74 78 82 86 90 94
     4÷π = 1.27283479368898552763028955877501991659132774562422849516943825241995907438793488819711543252100078   α terms=  1 2 2  2  2  2   2   2   2   2   2   2   2   2   2   2   2    2    2
                                                                                                                  ß terms=    1 9 25 49 81 121 169 225 289 361 441 529 625 729 841 961 1089 1225
 ½π, ½pi = 1.57071779679495409624134340842172803914665374867756685353559415360226497865248110814032818773478787   α terms=  1 ½ 1/3 ¼ 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17
   π, pi = 3.14259165433954305090112773725220456615353825631695587367530386050342717161955770321769607013860474   α terms=  0 1 2 2  2  2  2   2   2   2   2   2   2   2   2   2   2   2    2    2
                                                                                                                  ß terms=    4 1 9 25 49 81 121 169 225 289 361 441 529 625 729 841 961 1089 1225
   π, pi = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706   α terms=  0 1 3 5 7  9 11 13 15 17 19  21  23  25  27  29  31  33  35  37  39  41
                                                                                                                  ß terms=    4 1 4 9 16 25 36 49 64 81 100 121 144 169 196 225 256 289 324 361 400
   π, pi = 3.14159265358979298847014326453044038404101783047277203674633230347271153796007366409681897722403708   α terms=  3 6 6  6  6  6   6   6   6   6   6   6   6   6   6   6   6    6    6
                                                                                                                  ß terms=    1 9 25 49 81 121 169 225 289 361 441 529 625 729 841 961 1089 1225

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

Version 2 derived from PL/I

<lang rexx>/* REXX **************************************************************

Derived from PL/I with a little "massage"
SQRT2= 1.41421356237309505 <- PL/I Result
1.41421356237309504880168872421 <- REXX Result 30 digits
NAPIER= 2.71828182845904524
2.71828182845904523536028747135
PI= 3.14159262280484695
3.14159262280484694855146925223
06.09.2012 Walter Pachl
- - - - /

 Numeric Digits 30
 Parse Value '1 2 3 0 0' with Sqrt2 napier pi a b
 Say left('SQRT2=' ,10) calc(sqrt2,  200)
 Say left('NAPIER=',10) calc(napier, 200)
 Say left('PI='    ,10) calc(pi,     200)
 Exit

Get_Coeffs: procedure Expose a b Sqrt2 napier pi

 Parse Arg form, n
 select
   when form=Sqrt2 Then do
     if n > 0 then a = 2; else a = 1
     b = 1
     end
   when form=Napier Then do
     if n > 0 then a = n; else a = 2
     if n > 1 then b = n - 1; else b = 1
     end
   when form=pi Then do
     if n > 0 then a = 6; else a = 3
     b = (2*n - 1)**2
     end
   end
 Return

Calc: procedure Expose a b Sqrt2 napier pi

 Parse Arg form,n
 Temp=0
 do ni = n to 1 by -1
   Call Get_Coeffs form, ni
   Temp = B/(A + Temp)
   end
 call Get_Coeffs  form, 0
 return (A + Temp)</lang>

Version 3 better approximation

<lang rexx>/* REXX *************************************************************

The task description specifies a continued fraction for pi
that gives a reasonable approximation.
Literature shows a better CF that yields pi with a precision of
200 digits.
http://de.wikipedia.org/wiki/Kreiszahl
1
pi = 3 + ------------------------
1
7 + --------------------
1
15 + ---------------
1
1 + -----------
292 + ...
This program uses that CF and shows the first 50 digits
PI =3.1415926535897932384626433832795028841971693993751...
PIX=3.1415926535897932384626433832795028841971693993751...
201 correct digits
18.09.2012 Walter Pachl
- - - - /

 pi='3.1415926535897932384626433832795028841971'||,
    '693993751058209749445923078164062862089986280348'||,
    '253421170679821480865132823066470938446095505822'||,
    '317253594081284811174502841027019385211055596446'||,
    '229489549303819644288109756659334461284756482337'||,
    '867831652712019091456485669234603486104543266482'||,
    '133936072602491412737245870066063155881748815209'||,
    '209628292540917153643678925903600113305305488204'||,
    '665213841469519415116094330572703657595919530921'||,
    '861173819326117931051185480744623799627495673518'||,
    '857527248912279381830119491298336733624'
 Numeric Digits 1000
 al='7 15 1 292 1 1 1 2 1 3 1 14 2 1 1 2 2 2 2 1 84 2',
    '1 1 15 3 13 1 4 2 6 6 99 1 2 2 6 3 5 1 1 6 8 1 7 1 2',
    '3 7 1 2 1 1 12 1 1 1 3 1 1 8 1 1 2 1 6 1 1 5 2 2 3 1',
    '2 4 4 16 1 161 45 1 22 1 2 2 1 4 1 2 24 1 2 1 3 1 2',
    '1 1 10 2 5 4 1 2 2 8 1 5 2 2 26 1 4 1 1 8 2 42 2 1 7',
    '3 3 1 1 7 2 4 9 7 2 3 1 57 1 18 1 9 19 1 2 18 1 3 7',
    '30 1 1 1 3 3 3 1 2 8 1 1 2 1 15 1 2 13 1 2 1 4 1 12',
    '1 1 3 3 28 1 10 3 2 20 1 1 1 1 4 1 1 1 5 3 2 1 6 1 4'
 a.=3
 Do i=1 By 1 while al<>
   Parse Var al a.i al
   End
 pix=calc(194)
 Do e=1 To length(pi)
   If substr(pix,e,1)<>substr(pi,e,1) Then Leave
   End
 Numeric Digits 50
 Say 'PI ='||(pi+0)||'...'
 Say 'PIX='||(pix+0)||'...'
 Say (e-1) 'correct digits'
 Exit

Get_Coeffs: procedure Expose a b a.

 Parse Arg n
 a=a.n
 b=1
 Return

Calc: procedure Expose a b a.

 Parse Arg n
 Temp=0
 do ni = n to 1 by -1
   Call Get_Coeffs ni
   Temp = B/(A + Temp)
   end
 call Get_Coeffs 0
 return (A + Temp)</lang>

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