Continued fraction

(Redirected from Continued fractions)
Continued fraction
You are encouraged to solve this task according to the task description, using any language you may know.
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 a0 = 1 then aN = 2. bN is always 1.

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

For Napier's Constant, use a0 = 2, then aN = N. b1 = 1 then bN = N − 1.

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

For Pi, use a0 = 3 then aN = 6. bN = (2N − 1)2.

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

Contents

(The source text for these examples can also be found on Bitbucket.)

Generic function for estimating continued fractions:

generic   type Scalar is digits <>;    with function A (N : in Natural)  return Natural;   with function B (N : in Positive) return Natural;function Continued_Fraction (Steps : in Natural) return Scalar;
function Continued_Fraction (Steps : in Natural) return Scalar is   function A (N : in Natural)  return Scalar is (Scalar (Natural'(A (N))));   function B (N : in Positive) return Scalar is (Scalar (Natural'(B (N))));    Fraction : Scalar := 0.0;begin   for N in reverse Natural range 1 .. Steps loop      Fraction := B (N) / (A (N) + Fraction);   end loop;   return A (0) + Fraction;end Continued_Fraction;

Test program using the function above to estimate the square root of 2, Napiers constant and pi:

with Ada.Text_IO; with Continued_Fraction; procedure Test_Continued_Fractions is   type Scalar is digits 15;    package Square_Root_Of_2 is      function A (N : in Natural)  return Natural is (if N = 0 then 1 else 2);      function B (N : in Positive) return Natural is (1);       function Estimate is new Continued_Fraction (Scalar, A, B);   end Square_Root_Of_2;    package Napiers_Constant is      function A (N : in Natural)  return Natural is (if N = 0 then 2 else N);      function B (N : in Positive) return Natural is (if N = 1 then 1 else N-1);       function Estimate is new Continued_Fraction (Scalar, A, B);   end Napiers_Constant;    package Pi is      function A (N : in Natural)  return Natural is  (if N = 0 then 3 else 6);      function B (N : in Positive) return Natural is ((2 * N - 1) ** 2);       function Estimate is new Continued_Fraction (Scalar, A, B);   end Pi;    package Scalar_Text_IO is new Ada.Text_IO.Float_IO (Scalar);   use Ada.Text_IO, Scalar_Text_IO;begin   Put (Square_Root_Of_2.Estimate (200), Exp => 0); New_Line;   Put (Napiers_Constant.Estimate (200), Exp => 0); New_Line;   Put (Pi.Estimate (10000),             Exp => 0); New_Line;end Test_Continued_Fractions;

This example is exactly the same as the preceding one, but implemented using only Ada 95 features.

generic   type Scalar is digits <>;    with function A (N : in Natural)  return Natural;   with function B (N : in Positive) return Natural;function Continued_Fraction_Ada95 (Steps : in Natural) return Scalar;
function Continued_Fraction_Ada95 (Steps : in Natural) return Scalar is   function A (N : in Natural)  return Scalar is   begin      return Scalar (Natural'(A (N)));   end A;    function B (N : in Positive) return Scalar is   begin      return Scalar (Natural'(B (N)));   end B;    Fraction : Scalar := 0.0;begin   for N in reverse Natural range 1 .. Steps loop      Fraction := B (N) / (A (N) + Fraction);   end loop;   return A (0) + Fraction;end Continued_Fraction_Ada95;
with Ada.Text_IO; with Continued_Fraction_Ada95; procedure Test_Continued_Fractions_Ada95 is   type Scalar is digits 15;    package Square_Root_Of_2 is      function A (N : in Natural)  return Natural;      function B (N : in Positive) return Natural;       function Estimate is new Continued_Fraction_Ada95 (Scalar, A, B);   end Square_Root_Of_2;    package body Square_Root_Of_2 is      function A (N : in Natural) return Natural is      begin         if N = 0 then            return 1;         else            return 2;         end if;      end A;       function B (N : in Positive) return Natural is      begin         return 1;      end B;   end Square_Root_Of_2;    package Napiers_Constant is      function A (N : in Natural)  return Natural;      function B (N : in Positive) return Natural;       function Estimate is new Continued_Fraction_Ada95 (Scalar, A, B);   end Napiers_Constant;    package body Napiers_Constant is      function A (N : in Natural) return Natural is      begin         if N = 0 then            return 2;         else            return N;         end if;      end A;       function B (N : in Positive) return Natural is      begin          if N = 1 then             return 1;          else             return N - 1;          end if;      end B;   end Napiers_Constant;    package Pi is      function A (N : in Natural)  return Natural;      function B (N : in Positive) return Natural;       function Estimate is new Continued_Fraction_Ada95 (Scalar, A, B);   end Pi;    package body Pi is      function A (N : in Natural) return Natural is      begin         if N = 0 then            return 3;         else            return 6;         end if;      end A;       function B (N : in Positive) return Natural is      begin         return (2 * N - 1) ** 2;      end B;   end Pi;    package Scalar_Text_IO is new Ada.Text_IO.Float_IO (Scalar);   use Ada.Text_IO, Scalar_Text_IO;begin   Put (Square_Root_Of_2.Estimate (200), Exp => 0); New_Line;   Put (Napiers_Constant.Estimate (200), Exp => 0); New_Line;   Put (Pi.Estimate (10000),             Exp => 0); New_Line;end Test_Continued_Fractions_Ada95;
Output:
 1.41421356237310
2.71828182845905
3.14159265358954

ATS

A fairly direct translation of the C version without using advanced features of the type system:

(* a coefficient function creates double values from in paramters *)typedef coeff_f = int -> double (* a continued fraction is described by a record of two coefficent functions a and b *)typedef frac = @{a= coeff_f, b= coeff_f} (* recursive definition of the approximation *)fun calc_rec(f: frac, n: int, r: double): double =  if n = 0 then    f.a(0) + r // base case  else let // recursive case    val a: double = f.a(n) // record access    val b: double = f.b(n)    val s: double = b / (a + r)    //val () = printf("r%d = %f, a = %f, b = %f, r%d = %f\n", @(n, r, a, b, (n-1), s))  in calc_rec(f, n - 1, s) end fun calc(f: frac, n: int): double =  calc_rec(f, n, 0.0) (* anonymous coefficient functions created with lam (short for lambda) *)val sqrt2 = @{ a= lam (n: int): double => if n = 0 then 1.0 else 2.0,               b= lam (n: int): double => 1.0 } val napier = @{ a= lam (n: int): double => if n = 0 then 2.0 else 1.0 * (n),                b= lam (n: int): double => if n = 1 then 1.0 else n - 1.0 } fun square(x: double): double = x * x val pi = @{ a= lam (n: int): double => if n = 0 then 3.0 else 6.0,            b= lam (n: int): double => square (2.0 * n - 1) } implement main () = begin  printf ("sqrt2  = %f\n", @(calc(sqrt2,  100)));  printf ("napier = %f\n", @(calc(napier, 100)));  printf ("  pi   = %f\n", @(calc(  pi  , 100)));end

Axiom

Axiom provides a ContinuedFraction domain:

get(obj) == convergents(obj).1000 -- utility to extract the 1000th valueget continuedFraction(1, repeating [1], repeating [2]) :: Floatget continuedFraction(2, cons(1,[i for i in 1..]), [i for i in 1..]) :: Floatget continuedFraction(3, [i^2 for i in 1.. by 2], repeating [6]) :: Float
Output:
   (1)  1.4142135623 730950488                                                                  Type: Float    (2)  2.7182818284 590452354                                                                  Type: Float    (3)  3.1415926538 39792926                                                                  Type: Float

The value for π 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:

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+tempcf(1, repeating [1], repeating [2], 1000) :: Floatcf(2, cons(1,[i for i in 1..]), [i for i in 1..], 1000) :: Floatcf(3, [i^2 for i in 1.. by 2], repeating [6], 1000) :: Float

   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, e1.41421356237309504880168872420969807856967187537694807317667973799073247846212055511094575957753221653.1415924109 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274 Mathematica sqrt2=Function[n,{1,Transpose@{Array[2&,n],Array[1&,n]}}];napier=Function[n,{2,Transpose@{Range[n],Prepend[Range[n-1],1]}}];pi=Function[n,{3,Transpose@{Array[6&,n],Array[(2#-1)^2&,n]}}];approx=Function[l, N[Divide@@First@Fold[{{#2.#[[;;,1]],#2.#[[;;,2]]},#[[1]]}&,{{l[[2,1,1]]l[[1]]+l[[2,1,2]],l[[2,1,1]]},{l[[1]],1}},l[[2,2;;]]],10]];r2=approx/@{sqrt2@#,napier@#,pi@#}&@10000;r2//TableForm Output: 1.414213562 2.718281828 3.141592654  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 */ 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) 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 approximation. One should, therefore, not SHOW more that 15 digits. See http://de.wikipedia.org/wiki/Kreiszahl See Rexx for a better computation OCaml let pi = 3, fun n -> ((2*n-1)*(2*n-1), 6)and nap = 2, fun n -> (max 1 (n-1), n)and root2 = 1, fun n -> (1, 2) in let eval (i,f) k = let rec frac n = let a, b = f n in float a /. (float b +. if n >= k then 0.0 else frac (n+1)) in float i +. frac 1 in Printf.printf "sqrt(2)\t= %.15f\n" (eval root2 1000);Printf.printf "e\t= %.15f\n" (eval nap 1000);Printf.printf "pi\t= %.15f\n" (eval pi 1000); Output (inaccurate due to too few terms): sqrt(2) = 1.414213562373095 e = 2.718281828459046 pi = 3.141592653340542 PARI/GP Partial solution for simple continued fractions. back(v)=my(t=contfracpnqn(v));t[1,1]/t[2,1]*1.back(vector(100,i,2-(i==1))) Output: %1 = 1.4142135623730950488016887242096980786 Perl We'll use closures to implement the infinite lists of coeffficients. sub continued_fraction { my ($a, $b,$n) = (@_[0,1], $_[2] // 100);$a->() + ($n &&$b->() / continued_fraction($a,$b, $n-1));} printf "√2 ≈ %.9f\n", continued_fraction do { my$n; sub { $n++ ? 2 : 1 } }, sub { 1 };printf "e ≈ %.9f\n", continued_fraction do { my$n; sub { $n++ || 2 } }, do { my$n; sub { $n++ || 1 } };printf "π ≈ %.9f\n", continued_fraction do { my$n; sub { $n++ ? 6 : 3 } }, do { my$n; sub { (2*$n++ + 1)**2 } }, 1_000;printf "π/2 ≈ %.9f\n", continued_fraction do { my$n; sub { 1/($n++ || 1) } }, sub { 1 }, 1_000; Output: √2 ≈ 1.414213562 e ≈ 2.718281828 π ≈ 3.141592653 π/2 ≈ 1.570717797 Perl 6 sub continued-fraction(:@a, :@b, Int :$n = 100){    my $x = @a[$n - 1];    $x = @a[$_ - 1] + @b[$_] /$x for reverse 1 ..^ $n;$x;} printf "√2 ≈ %.9f\n", continued-fraction(:a(1, 2 xx *), :b(*, 1 xx *));printf "e  ≈ %.9f\n", continued-fraction(:a(2, 1 .. *), :b(*, 1, 1 .. *));printf "π  ≈ %.9f\n", continued-fraction(:a(3, 6 xx *), :b(*, [\+] 1, (8, 16 ... *)), :n(1000));
Output:
√2 ≈ 1.414213562
e  ≈ 2.718281828
π  ≈ 3.141592654

A more original and a bit more abstract method would consist in viewing a continued fraction on rank n as a function of a variable x:

$\mathrm{CF}_3(x) = a_0 + \cfrac{b_1}{a_1 + \cfrac{b_2}{a_2 + \cfrac{b_3}{a_3 + x}}}$

Or, more consistently:

$\mathrm{CF}_3(x) = a_0 + \cfrac{b_0}{a_1 + \cfrac{b_1}{a_2 + \cfrac{b_2}{a_3 + \cfrac{b_3}{x}}}}$

Viewed as such, CFn(x) could be written recursively:

$\mathrm{CF}_n(x) = \mathrm{CF}_{n-1}(a_n + \frac{b_n}{x})$

Or in other words:

$\mathrm{CF}_n= \mathrm{CF}_{n-1}\circ f_n = \mathrm{CF}_{n-2}\circ f_{n-1}\circ f_n=\ldots=f_0\circ f_1 \ldots \circ f_n$

where $f_n(x) = a_n + \frac{b_n}{x}$

Perl6 allows us to define a custom composition operator. We can then use it with the triangular reduction metaoperator, and map each resulting function with an infinite value for x (any value would do actually, but infinite make it consistent with this particular task).

sub infix:<⚬>(&f, &g) { -> $x { &f(&g($x)) } }sub continued-fraction(@a, @b) {    [\⚬] map { @a[$_] + @b[$_] / * }, 0 .. *} printf "√2 ≈ %.9f\n", continued-fraction((1, 2 xx *), (1 xx *))[10](Inf);printf "e  ≈ %.9f\n", continued-fraction((2, 1 .. *), (1, 1 .. *))[10](Inf);printf "π  ≈ %.9f\n", continued-fraction((3, 6 xx *), ((1, 3, 5 ... *) X** 2))[1000](Inf);

The main advantage is that the definition of the function does not need to know for which rank n it is computed. This is arguably closer to the mathematical definition.

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;
Output:
 1.41421356237309505E+0000 

Version for NAPIER:

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;
 2.71828182845904524E+0000 

Version for SQRT2, NAPIER, PI

/* 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;
Output:
SQRT2=     1.41421356237309505
NAPIER=    2.71828182845904524
PI=        3.14159265358979349


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 fractionscontinued_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 2sqrt_2_ab(0, 1, 1).sqrt_2_ab(_, 2, 1). % definitions for napiernapier_ab(0, 2, _).napier_ab(1, 1, 1).napier_ab(N, N, V) :-	V is N - 1. % definitions for pipi_ab(0, 3, _).pi_ab(N, 6, V) :-	V is (2 * N - 1)*(2 * N - 1).
Output:
 ?- continued_fraction.
sqrt(2) = 1.4142135623730951
e       = 2.7182818284590455
pi      = 3.141592622804847
true .


Python

Works with: Python version 2.6+ and 3.x
from fractions import Fractionimport itertoolstry: zip = itertools.izipexcept: pass # The Continued Fractiondef 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 stringdef 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 sqrt2def 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 Constantdef 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 Pidef 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

Fast iterative version

Translation of: D
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 = 50print calc(fsqrt2, 200)print calc(fnapier, 200)print calc(fpi, 200)
Output:
1.4142135623730950488016887242096980785696718753770
2.7182818284590452353602874713526624977572470937000
3.1415926839198062649342019294083175420335002640134

Racket

Using Doubles

This version uses standard double precision floating point numbers:

 #lang racket(define (calc cf n)  (match/values (cf 0)    [(a0 b0)     (+ a0        (for/fold ([t 0.0]) ([i (in-range (+ n 1) 0 -1)])          (match/values (cf i)                        [(a b) (/ b (+ a t))])))])) (define (cf-sqrt i)   (values  (if (> i 0) 2 1)  1))(define (cf-napier i) (values  (if (> i 0) i 2)  (if (> i 1) (- i 1) 1)))(define (cf-pi i)     (values  (if (> i 0) 6 3)  (sqr (- (* 2 i) 1)))) (calc cf-sqrt   200)(calc cf-napier 200)(calc cf-pi     200) 

Output:

 1.41421356237309512.71828182845904553.1415926839198063 

Version - Using Doubles

This versions uses big floats (arbitrary precision floating point):

 #lang racket(require math)(bf-precision 2048) ; in bits (define (calc cf n)  (match/values (cf 0)    [(a0 b0)     (bf+ (bf a0)        (for/fold ([t (bf 0)]) ([i (in-range (+ n 1) 0 -1)])          (match/values (cf i)                        [(a b) (bf/ (bf b) (bf+ (bf a) t))])))])) (define (cf-sqrt i)   (values  (if (> i 0) 2 1)  1))(define (cf-napier i) (values  (if (> i 0) i 2)  (if (> i 1) (- i 1) 1)))(define (cf-pi i)     (values  (if (> i 0) 6 3)  (sqr (- (* 2 i) 1)))) (calc cf-sqrt   200)(calc cf-napier 200)(calc cf-pi     200) 

Output:

 (bf #e1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358960036439214262599769155193770031712304888324413327207659690547583107739957489062466508437105234564161085482146113860092820802430986649987683947729823677905101453725898480737256099166805538057375451207262441039818826744940289448489312217214883459060818483750848688583833366310472320771259749181255428309841375829513581694269249380272698662595131575038315461736928338289219865139248048189188905788104310928762952913687232022557677738108337499350045588767581063729)(bf #e2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526059563073813232862794349076323382988075319525101901157383418793070215408914993488416750924476146066808226480016847741185374234544243710753907774499206955170276183860626133138458300075204493382656029760673711320070932870912744374704723624212700454495421842219077173525899689811474120614457405772696521446961165559468253835854362096088934714907384964847142748311021268578658461064714894910680584249490719358138073078291397044213736982988247857479512745588762993966446075)(bf #e3.14159268391980626493420192940831754203350026401337226640663040854412059241988978103217808449508253393479795573626200366332733859609651462659489470805432281782785922056335606047700127154963266242144951481397480765182268219697420028007903565511884267297358842935537138583640066772149177226656227031792115896439889412205871076985598822285367358003457939603015797225018209619662200081521930463480571130673429337524564941105654923909951299948539893933654293161126559643573974163405197696633200469475250152247413175932572922175467223988860975105100904322239324381097207835036465269418118204894206705789759765527734394105147) 

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.

/*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;   endN=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=10000a=1 rep(2)      do j=1 for T by 2;   b=b j**2;  end;              call tell '4÷π'/*══════════════════════════════════════════════════════════════════════*/T=10000a=1;  do j=1 for T;        a=a 1/j; @=@ '1/'j; end;     call tell '½π, ½pi'/*══════════════════════════════════════════════════════════════════════*/T=10000a=0 1 rep(2)      do j=1 for T by 2;   b=b j**2;     end;   b=4 b;  call tell 'π, pi'/*══════════════════════════════════════════════════════════════════════*/T=10000a=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=-.25          then w='-¼'         z=z w         endreturn 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

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

/* 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)

Version 3 better approximation

/* 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)

Ruby

require 'bigdecimal' # square root of 2sqrt2 = Object.newdef sqrt2.a(n); n == 1 ? 1 : 2; enddef sqrt2.b(n); 1; end # Napier's constantnapier = Object.newdef napier.a(n); n == 1 ? 2 : n - 1; enddef napier.b(n); n == 1 ? 1 : n - 1; end pi = Object.newdef pi.a(n); n == 1 ? 3 : 6; enddef 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  endend puts estimate(sqrt2, 50).to_s('F')puts estimate(napier, 50).to_s('F')puts estimate(pi, 10).to_s('F')
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. 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() }} 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: 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() }} 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. package require Tcl 8.6 # Term generators; yield list of pairsproc 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 calculatorproc 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} # Demonstrationputs [cf r2]puts [cf e]puts [cf pi 250]; # Converges more slowly
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.

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);]
Output:
16
1.414213562372820
1.414213562373100
13
2.718281828459380
2.718281828459050
10000
3.141592653589540
3.141592653589790
`