Continued fraction

Continued fraction
You are encouraged to solve this task according to the task description, using any language you may know.

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

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

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

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

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

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

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

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

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

Ada

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

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

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

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

begin

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

end ContFract; </lang>

Output:

 1.41421356237310
 2.71828182845905
 3.14159262280485

D

<lang d>import std.typecons;

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

   FP temp = 0.0;

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

}

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

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

}

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

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

}

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

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

}

void main() {

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

}</lang>

Output:

1.4142135623730950487
2.7182818284590452354
3.1415926228048469486

Factor

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

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

 math.ranges prettyprint sequences ;

IN: rosetta.cfrac

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

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

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

M: sqrt2 cfrac-b

   ! Always b_n is 1.
   2drop 1 ;

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

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

M: napier cfrac-b

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

SINGLETON: pi M: pi cfrac-a

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

M: pi cfrac-b

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

cfrac-estimate ( cfrac terms -- number )

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

decimalize ( rational prec -- string )

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

<PRIVATE

main ( -- )

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

PRIVATE>

MAIN: main</lang>

Output:

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

Go

<lang go>package main

import "fmt"

type cfTerm struct {

   a, b int

}

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

func cfSqrt2(nTerms int) cf {

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

}

func cfNap(nTerms int) cf {

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

}

func cfPi(nTerms int) cf {

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

}

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

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

}

func main() {

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

}</lang>

Output:

sqrt2: 1.4142135623730965
nap:   2.7182818284590455
pi:    3.141623806667839

Haskell

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

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

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

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

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

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

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

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

Output:

1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572735013846230912297024924836055850737212644121497099935831413222665927505592755799950501152782060571
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019
3.141592653297590947683406834261190738869139611505752231394089152890909495973464508817163306557131591579057202097715021166512662872910519439747609829479577279606075707015622200744006783543589980682386

J

<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

  decimalize=: 4 : 0"0 NB. decimalize translated from python
'q y' =. y (<.@% , |~) 1
res=. ":q
res=. res,'.'
for. i.x do.
 y=. 10*y
 'q y'=. y (<.@% , |~) 1
 res=. res,":q
end.

)

  30 10 30 decimalize sqrt2, pi, e

1.414213562373095048801688724209 3.1415924109 2.718281828459045235360287471352 </lang>

Python

Works with: Python version 2.6+ and 3.x

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

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

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

1. 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. 1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273501384623091229702492483605585073721264412149709993583141322266592750559275579995050115278206057147

1. 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))

1. 2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526059563073813232862794349076323382988075319525101901

1. 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))

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

Ruby

<lang ruby>require 'bigdecimal'

1. square root of 2

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

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

1. Estimates the value of a continued fraction _cfrac_, to _prec_
2. decimal digits of precision. Returns a BigDecimal. _cfrac_ must
3. 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

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

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

   }

}

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

}

1. Demonstration

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

Output:

1.4142135623730951
2.7182818284590455
3.1415926373965735