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>

BBC BASIC

Works with: BBC BASIC for Windows

<lang bbcbasic> *FLOAT64

     @% = &1001010
     
     PRINT "SQR(2) = " ; FNcontfrac(1, 1, "2", "1")
     PRINT "     e = " ; FNcontfrac(2, 1, "N", "N")
     PRINT "    PI = " ; FNcontfrac(3, 1, "6", "(2*N+1)^2")
     END
     
     REM a$ and b$ are functions of N
     DEF FNcontfrac(a0, b1, a$, b$)
     LOCAL N, expr$
     REPEAT
       N += 1
       expr$ += STR$(EVAL(a$)) + "+" + STR$(EVAL(b$)) + "/("
     UNTIL LEN(expr$) > (65500 - N)
     = a0 + b1 / EVAL (expr$ + "1" + STRING$(N, ")"))</lang>

Output:

SQR(2) = 1.414213562373095
     e = 2.718281828459046
    PI = 3.141592653588017

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

Common Lisp

Translation of: C++

<lang lisp>(defun estimate-continued-fraction (generator n)

 (let ((temp 0))
   (loop for n1 from n downto 1
      do (multiple-value-bind (a b)

(funcall generator n1) (setf temp (/ b (+ a temp)))))

   (+ (funcall generator 0) temp)))

(format t "sqrt(2) = ~a~%" (coerce (estimate-continued-fraction (lambda (n) (values (if (> n 0) 2 1) 1)) 20) 'double-float)) (format t "napier's = ~a~%" (coerce (estimate-continued-fraction (lambda (n) (values (if (> n 0) n 2) (if (> n 1) (1- n) 1))) 15) 'double-float))

(format t "pi = ~a~%" (coerce (estimate-continued-fraction (lambda (n) (values (if (> n 0) 6 3) (* (1- (* 2 n)) (1- (* 2 n))))) 10000) 'double-float)) </lang>

Output:

sqrt(2) = 1.4142135623730947d0
napier's = 2.7182818284590464d0
pi = 3.141592653589543d0

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>

Mathematica

<lang 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}&,{{l2,1,1l 1+l2,1,2,l2,1,1},{l1,1}},l2,2;;],10]]; r2=approx/@{sqrt2@#,napier@#,pi@#}&@10000;r2//TableForm </lang> Output:

1.414213562
2.718281828
3.141592654

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

Perl

Translation of: Perl6

Perl5 does not really support lazy lists but we can use closures to do something similar.

With a recent version of Perl (5.10 or above), we can also use persistent private variables (see perlsub). We just need to declare this with the 'use feature "state"' pragma.

In this example we'll show both methods.

<lang perl>sub continued_fraction {

   my ($a, $b, $n) = @_;
   $n //= 100;

   my @a = map &$a, 0 .. $n;
   my @b = map &$b, 1 .. $n;

   my $x = pop @a;
   $x = pop(@a) + pop(@b)/$x while @a;
   return $x;

}

using closures

printf "√2 ≈ %.9f\n", continued_fraction do { my $n; sub { $n++ ? 2 : 1 } }, sub { 1 };

using persistent private variables

use feature 'state'; printf "e ≈ %.9f\n", continued_fraction sub { state $n; $n++ ? $n-1 : 2 }, sub { state $n; $n++ ? $n-1 : 1 }; printf "τ ≈ %.9f\n", continued_fraction sub { 6 }, sub { state $n; $n++ ? (2*$n - 1)**2 : 2 }, 1_000; </lang>

Output:

√2 ≈ 1.414213562
e  ≈ 2.718281828
τ  ≈ 6.283185307

Perl 6

<lang perl6>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));</lang>

Output:

√2 ≈ 1.414213562
e  ≈ 2.718281828
π  ≈ 3.141592654

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=-.25          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