# 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:
${\displaystyle 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 ${\displaystyle a_{0}=1}$ then ${\displaystyle a_{N}=2}$. ${\displaystyle b_{N}}$ is always ${\displaystyle 1}$.

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

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

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

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

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

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

### Using only Ada 95 features

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

## ALGOL 68

Works with: Algol 68 Genie version 2.8
 PROC cf = (INT steps, PROC (INT) INT a, PROC (INT) INT b) REAL:BEGIN  REAL result;  result := 0;  FOR n FROM steps BY -1 TO 1 DO      result := b(n) / (a(n) + result)  OD;  a(0) + resultEND; PROC asqr2 = (INT n) INT: (n = 0 | 1 | 2);PROC bsqr2 = (INT n) INT: 1; PROC anap = (INT n) INT: (n = 0 | 2 | n);PROC bnap = (INT n) INT: (n = 1 | 1 | n - 1); PROC api = (INT n) INT: (n = 0 | 3 | 6);PROC bpi = (INT n) INT: (n = 1 | 1 | (2 * n - 1) ** 2); INT precision = 10000; print (("Precision: ", precision, newline));print (("Sqr(2):    ", cf(precision, asqr2, bsqr2), newline));print (("Napier:    ", cf(precision, anap, bnap), newline));print (("Pi:        ", cf(precision, api, bpi)))
Output:

Precision:      +10000
Sqr(2):    +1.41421356237310e  +0
Napier:    +2.71828182845905e  +0
Pi:        +3.14159265358954e  +0


## ATS

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

#include"share/atspre_staload.hats"//(* ****** ****** *)//(*** 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}//(* ****** ****** *) fun calc(  f: frac, n: int) : double = let//(*** recursive definition of the approximation*)fun loop(  n: int, r: double) : double =(if n = 0  then f.a(0) + r  else loop (n - 1, f.b(n) / (f.a(n) + r))// end of [if])//in  loop (n, 0.0)end // end of [calc] (* ****** ****** *) val sqrt2 = @{  a= lam (n: int): double => if n = 0 then 1.0 else 2.0,  b= lam (n: int): double => 1.0} (* end of [val] *) 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} (* end of [val] *) val pi = @{  a= lam (n: int): double => if n = 0 then 3.0 else 6.0,  b= lam (n: int): double => let val x = 2.0 * n - 1 in x * x end} (* ****** ****** *) implementmain0 () =(  println! ("sqrt2  = ", calc(sqrt2,  100));  println! ("napier = ", calc(napier, 100));  println! ("  pi   = ", calc(  pi  , 100));) (* end of [main0] *)

## AutoHotkey

sqrt2_a(n) ; function definition is as simple as that{return n?2.0:1.0} sqrt2_b(n){return 1.0} napier_a(n){return n?n:2.0} napier_b(n){return n>1.0?n-1.0:1.0} pi_a(n){return n?6.0:3.0} pi_b(n){return (2.0*n-1.0)**2.0 ; exponentiation operator} calc(f,expansions){r:=0,i:=expansions; A nasty trick: the names of the two coefficient functions are generated dynamically; a dot surrounded by spaces means string concatenationf_a:=f . "_a",f_b:=f . "_b" while i>0 {; You can see two dynamic function calls herer:=%f_b%(i)/(%f_a%(i)+r)i--} return %f_a%(0)+r} Msgbox, % "sqrt 2 = " . calc("sqrt2", 1000) . "ne = " . calc("napier", 1000) . "npi = " . calc("pi", 1000)

Output with Autohotkey v1 (currently 1.1.16.05):

sqrt 2 = 1.414214e = 2.718282pi = 3.141593

Output with Autohotkey v2 (currently alpha 56):

sqrt 2 = 1.4142135623730951e = 2.7182818284590455pi = 3.1415926533405418

Note the far superiour accuracy of v2.

## 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 ${\displaystyle \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:

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

## Forth

Translation of: D
: 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;
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


## Fortran

module continued_fractions  implicit none   integer, parameter :: long = selected_real_kind(7,99)   type continued_fraction    integer                            :: a0, b1    procedure(series), pointer, nopass :: a, b  end type   interface    pure function series (n)      integer, intent(in) :: n      integer             :: series    end function  end interface contains   pure function define_cf (a0,a,b1,b) result(x)    integer, intent(in)           :: a0    procedure(series)             :: a    integer, intent(in), optional :: b1    procedure(series),   optional :: b    type(continued_fraction)      :: x    x%a0 = a0    x%a => a    if ( present(b1) ) then       x%b1 = b1    else       x%b1 = 1    end if    if ( present(b) ) then       x%b => b    else       x%b => const_1    end if  end function define_cf   pure integer function const_1(n)    integer,intent(in) :: n    const_1 = 1    end function   pure real(kind=long) function output(x,iterations)    type(continued_fraction), intent(in) :: x    integer,                  intent(in) :: iterations    integer                              :: i    output = x%a(iterations)    do i = iterations-1,1,-1      output = x%a(i) + (x%b(i+1) / output)    end do    output = x%a0 + (x%b1 / output)  end function output end module continued_fractions  program examples  use continued_fractions   type(continued_fraction) :: sqr2,napier,pi   sqr2   = define_cf(1,a_sqr2)  napier = define_cf(2,a_napier,1,b_napier)  pi     = define_cf(3,a=a_pi,b=b_pi)   write (*,*) output(sqr2,10000)  write (*,*) output(napier,10000)  write (*,*) output(pi,10000) contains   pure integer function a_sqr2(n)    integer,intent(in) :: n    a_sqr2 = 2  end function   pure integer function a_napier(n)    integer,intent(in) :: n    a_napier = n  end function   pure integer function b_napier(n)    integer,intent(in) :: n    b_napier = n-1  end function   pure integer function a_pi(n)    integer,intent(in) :: n    a_pi = 6  end function   pure integer function b_pi(n)    integer,intent(in) :: n    b_pi = (2*n-1)*(2*n-1)  end function end program examples
Output:
   1.4142135623730951
2.7182818284590455
3.1415926535895435


## 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())}
Output:
sqrt2: 1.4142135623730965
nap:   2.7182818284590455
pi:    3.141623806667839


import Data.List (unfoldr)import Data.Char (intToDigit) -- continued fraction represented as a (possibly infinite) list of pairssqrt2, 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 iterationsapproxCF :: (Integral a, Fractional b) => Int -> [(a, a)] -> bapproxCF t =  foldr (\(a,b) z -> fromIntegral a + fromIntegral b / z) 1 . take t -- infinite decimal representation of a real numberdecString :: RealFrac a => a -> StringdecString 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]
Output:
1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572735013846230912297024924836055850737212644121497099935831413222665927505592755799950501152782060571
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019
3.141592653297590947683406834261190738869139611505752231394089152890909495973464508817163306557131591579057202097715021166512662872910519439747609829479577279606075707015622200744006783543589980682386

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 termscf2rat n = foldr (\(a,b) f -> (a%1) + ((b%1) / f)) (1%1) . take n -- truncate after error is at most 1/pcf2rat_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 ## 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, e1.41421356237309504880168872420969807856967187537694807317667973799073247846212055511094575957753221653.1415924109                                                                                          2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274

Note that there are two kinds of continued fractions. The kind here where we alternate between a and b values, but in some other tasks b is always 1 (and not included in the list we use to represent the continued fraction). The other kind is evaluated in J using (+%)/ instead of +%/.

## jq

Works with: jq version 1.4

We take one of the points of interest here to be the task of representing the infinite series a0, a1, .... and b0, b1, .... compactly, preferably functionally. For the type of series typically encountered in continued fractions, this is most readily accomplished in jq 1.4 using a filter (a function), here called "next", which, given the triple [i, [a[i], b[i]], will produce the next triple [i+1, a[i+1], b[i+1]].

Another point of interest is avoiding having to specify the number of iterations. The approach adopted here allows one to specify the desired accuracy; in some cases, it is feasible to specify that the computation should continue until the accuracy permitted by the underlying floating point representation is achieved. This is done by specifying delta as 0, as shown in the examples.

We therefore proceed in two steps: continued_fraction( first; next; count ) computes an approximation based on the first "count" terms; and then continued_fraction_delta(first; next; delta) computes the continued fraction until the difference in approximations is less than or equal to delta, which may be 0, as previously noted.

 # "first" is the first triple,# e.g. [1,a,b]; count specifies the number of terms to use.def continued_fraction( first; next; count ):  # input: [i, a, b]]  def cf:      if .[0] == count then 0      else next as $ab | .[1] + (.[2] / ($ab | cf))      end ;  first | cf; # "first" and "next" are as above;# if delta is 0 then continue until there is no detectable change.def continued_fraction_delta(first; next; delta):  def abs: if . < 0 then -. else . end;  def cf:    # state: [n, prev]    .[0] as $n | .[1] as$prev    |  continued_fraction(first; next; $n+1) as$this    | if $prev == null then [$n+1, $this] | cf elif delta <= 0 and ($prev == $this) then$this      elif (($prev -$this)|abs) <= delta then $this else [$n+1, $this] | cf end; [2,null] | cf;  Examples: The convergence for pi is slow so we select delta = 1e-12 in that case. "Value : Direct : Continued Fraction", "2|sqrt : \(2|sqrt) : \(continued_fraction_delta( [1,1,1]; [.[0]+1, 2, 1]; 0))", "1|exp : \(1|exp) : \(2 + (1 / (continued_fraction_delta( [1,1,1]; [.[0]+1, .[1]+1, .[2]+1]; 0))))", "pi : \(1|atan * 4) : \(continued_fraction_delta( [1,3,1]; .[0]+1 | [., 6, ((2*. - 1) | (.*.))]; 1e-12)) (1e-12)"  Output: $ jq -M -n -r -f Continued_fraction.jqValue  :        Direct      : Continued Fraction2|sqrt : 1.4142135623730951 : 1.41421356237309511|exp  : 2.718281828459045  : 2.7182818284590455pi     : 3.141592653589793  : 3.1415926535892935 (1e-12)

function sqrt2(a,n)if areturn n>0?2.0:1.0elsereturn 1.0endend function napier(a,n)if areturn n>0?Float64(n):2.0elsereturn n>1?n-1.0:1.0endend function _pi(a,n)if areturn n>0?6.0:3.0elsereturn (2.0*n-1.0)^2.0 # exponentiation operatorendend function calc(f,expansions)a=true;b=falser=0.0for i=expansions:-1:1r=f(b,i)/(f(a,i)+r)endreturn f(a,0)+rend print("""sqrt 2 = $(calc(sqrt2, 1000))e =$(calc(napier, 1000)) pi = $(calc(_pi, 1000))""") Output: sqrt 2 = 1.4142135623730951e = 2.7182818284590455 pi = 3.141592653340542 ## Mathematica / Wolfram Language 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 Works with: rakudo version 2015-10-31 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(Nil, |(1 xx *)));printf "e  ≈%.9f\n", continued-fraction(:a(2, |(1 .. *)), :b(Nil, 1, |(1 .. *)));printf "π  ≈%.9f\n", continued-fraction(:a(3, |(6 xx *)), :b(Nil, |((1, 3, 5 ... *) X** 2)));
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:

${\displaystyle \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:

${\displaystyle \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, ${\displaystyle \mathrm {CF} _{n}(x)}$ could be written recursively:

${\displaystyle \mathrm {CF} _{n}(x)=\mathrm {CF} _{n-1}(a_{n}+{\frac {b_{n}}{x}})}$

Or in other words:

${\displaystyle \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 ${\displaystyle f_{n}(x)=a_{n}+{\frac {b_{n}}{x}}}$

Perl6 has a builtin composition operator. We can use it with the triangular reduction metaoperator, and evaluate each resulting function at infinity (any value would do actually, but infinite makes it consistent with this particular task).

sub continued-fraction(@a, @b) {    map { .(Inf) }, [\o] map { @a[$_] + @b[$_] / * }, ^Inf} printf "√2 ≈ %.9f\n", continued-fraction((1, |(2 xx *)), (1 xx *))[10];printf "e  ≈ %.9f\n", continued-fraction((2, |(1 .. *)), (1, |(1 .. *)))[10];printf "π  ≈ %.9f\n", continued-fraction((3, |(6 xx *)), ((1, 3, 5 ... *) X** 2))[100];
Output:
√2 ≈ 1.414213552
e  ≈ 2.718281827
π  ≈ 3.141592411

## 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 ## Sidef Translation of: Perl func continued_fraction(a, _, (0)) { a() }func continued_fraction(a, b, n=100) { a() + (b() / continued_fraction(a, b, n-1));} var f = Hash.new( "√2" => [ do { var n = 0; { n++ ? 2 : 1 } }, { 1 }, ], "e" => [ do { var n = 0; { n++ || 2 } }, do { var n = 0; { n++ || 1 } }, ], "π" => [ do { var n = 0; { n++ ? 6 : 3 } }, do { var n = 0; { (2*(n++) + 1)**2 } }, 1_000, ], "π/2" => [ do { var n = 0; { 1/(n++ || 1) } }, { 1 }, 1_000, ]); f.each { |k,v| printf("%3s ≈ %.9f\n", k, continued_fraction(v...));} Output:  e ≈ 2.718281828 π ≈ 3.141592653 √2 ≈ 1.414213562 π/2 ≈ 1.570011909  ## 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


## zkl

fcn cf(fa,fb,a0){fcn(fa,fb,a0,n){   a0 + [n..1,-1].reduce(        'wrap(p,n){ fb(n)/(fa(n)+p) },0.0) }.fp(fa,fb,a0)}

cf creates a function that calculates the continued fraction from the bottom up. The new function takes a single parameter, n, which is used to calculate the nth term.

sqrt2:=cf((2.0).noop,(1.0).noop,1.0);   sqrt2(200) : "%.20e".fmt(_).println();nap:=cf((0.0).create,fcn(n){ (n==1) and 1.0 or (n-1).toFloat() },2.0);   println(nap(15) - (1.0).e);pi:=cf((6.0).noop,fcn(n){ n=2*n-1; (n*n).toFloat() },3.0);   println(pi(1000) - (1.0).pi);

(1.0).create(n) --> n, (1.0).noop(n) --> 1.0

Output:
1.41421356237309514547e+00
1.33227e-15
-2.49251e-10