Continued fraction
A number may be represented as a continued fraction (see Mathworld for more information) as follows:
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.
For the square root of 2, is always . is then is .
For Napier's Constant, is , then is . is then is .
For Pi, is 3 then is always . is .
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
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
Python
<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 while True: yield 2
def sqrt2_b():
while True: yield 1
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 n = 1 while True: yield n n += 1
def Napier_b():
n=1 yield n while True: yield n n += 1
cf = CF(Napier_a(), Napier_b(), 950) print(pRes(cf, 200))
- 2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526059563073813232862794349076323382988075319525101901
- Test the Continued Fraction for Pi
def Pi_a():
yield 3 while True: yield 6
def Pi_b():
x = 1 while True: yield x*x x += 2
cf = CF(Pi_a(), Pi_b(), 950) print(pRes(cf, 10))
- 3.1415926532</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| result = BigDecimal(cfrac.b(n)).div(result, prec) + cfrac.a(n) 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
Note that Tcl does not provide arbitrary precision floating point numbers by default, so all result computations are done with IEEE double
s.
<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