Continued fraction: Difference between revisions
(Simplify the re-implementation) |
|||
Line 65: | Line 65: | ||
=={{header|Axiom}}== |
=={{header|Axiom}}== |
||
Axiom provides a ContinuedFraction domain: |
Axiom provides a ContinuedFraction domain: |
||
<lang Axiom>get(obj) == convergents(obj).1000 |
<lang Axiom>get(obj) == convergents(obj).1000 -- utility to extract the 1000th value |
||
get continuedFraction(1, repeating [1], repeating [2]) |
get continuedFraction(1, repeating [1], repeating [2]) :: Float |
||
get continuedFraction(2, cons(1,[i for i in 1..]), [i for i in 1..]) |
get continuedFraction(2, cons(1,[i for i in 1..]), [i for i in 1..]) :: Float |
||
get continuedFraction(3, [(2*i-1)^2 for i in 1..], repeating [6])</lang> |
get continuedFraction(3, [(2*i-1)^2 for i in 1..], repeating [6]) :: Float</lang> |
||
Output:<lang Axiom> (1) 1.4142135623 730950488 |
Output:<lang Axiom> (1) 1.4142135623 730950488 |
||
Type: Float |
Type: Float |
||
Line 79: | Line 79: | ||
The value of pi has an accuracy to only 9 decimal places after 1000 iterations, with an accuracy to 12 decimal places after 10000 iterations. |
The value of 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 |
We could re-implement this, with the same output:<lang Axiom>cf(initial, a, b, n) == |
||
⚫ | |||
n:PositiveInteger):Fraction(Integer) == |
|||
temp := 0 |
|||
⚫ | |||
temp : Fraction(Integer) := 0 |
|||
for i in (n-1)..1 by -1 repeat |
for i in (n-1)..1 by -1 repeat |
||
temp := a.i/(b.i+temp) |
temp := a.i/(b.i+temp) |
||
initial+temp |
initial+temp |
||
cf(1, repeating [1], repeating [2], 1000) :: Float |
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(2, cons(1,[i for i in 1..]), [i for i in 1..], 1000) :: Float |
||
cf(3, [(2*i-1)^2 for i in 1..], repeating [6], 1000) :: Float</lang> |
cf(3, [(2*i-1)^2 for i in 1..], repeating [6], 1000) :: Float</lang> |
||
Revision as of 21:29, 21 March 2012
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:
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 then . is always .
For Napier's Constant, use , then . then .
For Pi, use then . .
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
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, [(2*i-1)^2 for i in 1..], 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 of 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, [(2*i-1)^2 for i in 1..], repeating [6], 1000) :: Float</lang>
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
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>
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 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
<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'
- 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
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