Continued fraction: Difference between revisions
Added Forth |
m →{{header|REXX}}: deleted a blank line, expanded ''output'' a wee bit. -- ~~~~ |
||
Line 548:
/*also: [e+1] / [e-1] */
exit
/*────────────────────────────────$CF subroutine────────────────────────*/
$cf: procedure; parse arg c x,y; !=0; numeric digits digits()*2
Line 560 ⟶ 559:
divzero: say; say '***error!***'; say 'division by zero.'; say; exit 13
tell:parse arg _,v;say right(_,8) '=' format(v) 'a terms='left(a,40);a=;b=;return</lang>
'''output'''
<pre style="height:
sqrt(2) = 1.41421356237309504880168872420969807856967187537694807317668 a terms= 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
e = 2.71828182845904523536028747135266249775724709369995957496697 a terms= 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
|
Revision as of 15:12, 15 May 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, [i^2 for i in 1.. by 2], repeating [6]) :: Float</lang> Output:<lang Axiom> (1) 1.4142135623 730950488
Type: Float
(2) 2.7182818284 590452354 Type: Float
(3) 3.1415926538 39792926 Type: Float</lang>
The value for has an accuracy to only 9 decimal places after 1000 iterations, with an accuracy to 12 decimal places after 10000 iterations.
We could re-implement this, with the same output:<lang Axiom>cf(initial, a, b, n) ==
n=1 => initial temp := 0 for i in (n-1)..1 by -1 repeat temp := a.i/(b.i+temp) initial+temp
cf(1, repeating [1], repeating [2], 1000) :: Float cf(2, cons(1,[i for i in 1..]), [i for i in 1..], 1000) :: Float cf(3, [i^2 for i in 1.. by 2], repeating [6], 1000) :: Float</lang>
CoffeeScript
<lang coffeescript># Compute a continuous fraction of the form
- a0 + b1 / (a1 + b2 / (a2 + b3 / ...
continuous_fraction = (f) ->
a = f.a b = f.b c = 1 for n in [100000..1] c = b(n) / (a(n) + c) a(0) + c
- A little helper.
p = (a, b) ->
console.log a console.log b console.log "---"
do ->
fsqrt2 = a: (n) -> if n is 0 then 1 else 2 b: (n) -> 1 p Math.sqrt(2), continuous_fraction(fsqrt2)
fnapier = a: (n) -> if n is 0 then 2 else n b: (n) -> if n is 1 then 1 else n - 1 p Math.E, continuous_fraction(fnapier)
fpi = a: (n) -> return 3 if n is 0 6 b: (n) -> x = 2*n - 1 x * x p Math.PI, continuous_fraction(fpi)
</lang> output
> coffee continued_fraction.coffee 1.4142135623730951 1.4142135623730951 --- 2.718281828459045 2.7182818284590455 --- 3.141592653589793 3.141592653589793 ---
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
Forth
<lang forth>: fsqrt2 1 s>f 0> if 2 s>f else fdup then ;
- fnapier dup dup 1 > if 1- else drop 1 then s>f dup 1 < if drop 2 then s>f ;
- fpi dup 2* 1- dup * s>f 0> if 6 else 3 then s>f ;
( n -- f1 f2)
- cont.fraction ( xt n -- f)
1 swap 1+ 0 s>f \ calculate for 1 .. n do i over execute frot f+ f/ -1 +loop 0 swap execute fnip f+ \ calcucate for 0
- </lang>
- Output:
' fsqrt2 200 cont.fraction f. cr 1.4142135623731 ok ' fnapier 200 cont.fraction f. cr 2.71828182845905 ok ' fpi 200 cont.fraction f. cr 3.14159268391981 ok
Go
<lang go>package main
import "fmt"
type cfTerm struct {
a, b int
}
// follows subscript convention of mathworld and WP where there is no b(0). // cf[0].b is unused in this representation. type cf []cfTerm
func cfSqrt2(nTerms int) cf {
f := make(cf, nTerms) for n := range f { f[n] = cfTerm{2, 1} } f[0].a = 1 return f
}
func cfNap(nTerms int) cf {
f := make(cf, nTerms) for n := range f { f[n] = cfTerm{n, n - 1} } f[0].a = 2 f[1].b = 1 return f
}
func cfPi(nTerms int) cf {
f := make(cf, nTerms) for n := range f { g := 2*n - 1 f[n] = cfTerm{6, g * g} } f[0].a = 3 return f
}
func (f cf) real() (r float64) {
for n := len(f) - 1; n > 0; n-- { r = float64(f[n].b) / (float64(f[n].a) + r) } return r + float64(f[0].a)
}
func main() {
fmt.Println("sqrt2:", cfSqrt2(20).real()) fmt.Println("nap: ", cfNap(20).real()) fmt.Println("pi: ", cfPi(20).real())
}</lang>
- Output:
sqrt2: 1.4142135623730965 nap: 2.7182818284590455 pi: 3.141623806667839
Haskell
<lang haskell>import Data.List (unfoldr) import Data.Char (intToDigit)
-- continued fraction represented as a (possibly infinite) list of pairs sqrt2, napier, myPi :: [(Integer, Integer)] sqrt2 = zip (1 : [2,2..]) [1,1..] napier = zip (2 : [1..]) (1 : [1..]) myPi = zip (3 : [6,6..]) (map (^2) [1,3..])
-- approximate a continued fraction after certain number of iterations approxCF :: (Integral a, Fractional b) => Int -> [(a, a)] -> b approxCF t =
foldr (\(a,b) z -> fromIntegral a + fromIntegral b / z) 1 . take t
-- infinite decimal representation of a real number decString :: RealFrac a => a -> String decString frac = show i ++ '.' : decString' f where
(i,f) = properFraction frac decString' = map intToDigit . unfoldr (Just . properFraction . (10*))
main :: IO () main = mapM_ (putStrLn . take 200 . decString .
(approxCF 950 :: [(Integer, Integer)] -> Rational)) [sqrt2, napier, myPi]</lang>
- Output:
1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572735013846230912297024924836055850737212644121497099935831413222665927505592755799950501152782060571 2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019 3.141592653297590947683406834261190738869139611505752231394089152890909495973464508817163306557131591579057202097715021166512662872910519439747609829479577279606075707015622200744006783543589980682386
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>
Maxima
<lang maxima>cfeval(x) := block([a, b, n, z], a: x[1], b: x[2], n: length(a), z: 0,
for i from n step -1 thru 2 do z: b[i]/(a[i] + z), a[1] + z)$
cf_sqrt2(n) := [cons(1, makelist(2, i, 2, n)), cons(0, makelist(1, i, 2, n))]$
cf_e(n) := [cons(2, makelist(i, i, 1, n - 1)), append([0, 1], makelist(i, i, 1, n - 2))]$
cf_pi(n) := [cons(3, makelist(6, i, 2, n)), cons(0, makelist((2*i - 1)^2, i, 1, n - 1))]$
cfeval(cf_sqrt2(20)), numer; /* 1.414213562373097 */ % - sqrt(2), numer; /* 1.3322676295501878*10^-15 */
cfeval(cf_e(20)), numer; /* 2.718281828459046 */ % - %e, numer; /* 4.4408920985006262*10^-16 */
cfeval(cf_pi(20)), numer; /* 3.141623806667839 */ % - %pi, numer; /* 3.115307804568701*10^-5 */
/* convergence is much slower for pi */
fpprec: 20$
x: cfeval(cf_pi(10000))$
bfloat(x - %pi); /* 2.4999999900104930006b-13 */</lang>
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
REXX
The $CF subroutine (for continous fractions) isn't limited to positive integers,
any form of REXX numbers (negative, exponentiated, decimal fractions) can be used.
There isn't any realistic limit on the precision that can be used, although 100,000 digits would be a bit unwieldly to display.
<lang rexx>/*REXX program to calculate some contined fractions. */
numeric digits 60 /*let's use sixty digits for show&tell*/ terms=100 /* ... and an even one hundred terms. */ /*──────────────────────────────────────────────────────────────────────*/ a=copies(' 2',terms); call tell 'sqrt(2)',$cf(1 a) /*──────────────────────────────────────────────────────────────────────*/
do j=1 for terms a=a j end; call tell 'e',$cf(2 a,1 a)
/*──────────────────────────────────────────────────────────────────────*/ a=copies(' 6',terms)
do j=1 for terms b=b (2*j-1)**2 end; call tell 'pi',$cf(3 a,b)
/*──────────────────────────────────────────────────────────────────────*/ a=copies(' 1',terms); call tell 'phi',$cf(1 a) /*──────────────────────────────────────────────────────────────────────*/ a=copies(' 1 2',terms%2); call tell 'sqrt(3)',$cf(1 a)
/*also: 2∙sin(π/3) */
/*──────────────────────────────────────────────────────────────────────*/
do j=1 for terms%2 by 2 a=a j 1 end; call tell 'tan(1)',$cf(1 a)
/*──────────────────────────────────────────────────────────────────────*/
do j=1 for terms a=a 2*j+1 end; call tell 'coth(1)',$cf(1 a)
/*──────────────────────────────────────────────────────────────────────*/
do j=1 for terms a=a 4*j+2 end; call tell 'coth(½)',$cf(2 a) /*also: [e+1] / [e-1] */
exit /*────────────────────────────────$CF subroutine────────────────────────*/ $cf: procedure; parse arg c x,y; !=0; numeric digits digits()*2
do j=words(x) to 1 by -1 a=word(x,j); b=word(y,j); if b== then b=1 d=a+!; if d=0 then call divzero /*just in case divisor is bogus.*/ !=b/d end
return !+c
divzero: say; say '***error!***'; say 'division by zero.'; say; exit 13 tell:parse arg _,v;say right(_,8) '=' format(v) 'a terms='left(a,40);a=;b=;return</lang> output
sqrt(2) = 1.41421356237309504880168872420969807856967187537694807317668 a terms= 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 e = 2.71828182845904523536028747135266249775724709369995957496697 a terms= 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 pi = 3.14159241097198067426258886021672643729365090674537669567319 a terms= 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 phi = 1.61803398874989484820458683436563811772030781841853559947667 a terms= 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 sqrt(3) = 1.73205080756887729352744634150587236694280525381038062805521 a terms= 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 tan(1) = 1.55740772465490223050697480745836017308725077238152003838395 a terms= 1 1 3 1 5 1 7 1 9 1 11 1 13 1 15 1 17 1 coth(1) = 1.31303528549933130363616124693084783291201394124045265554315 a terms= 3 5 7 9 11 13 15 17 19 21 23 25 27 29 3 coth(½) = 2.16395341373865284877000401021802311709373860215079227253357 a terms= 6 10 14 18 22 26 30 34 38 42 46 50 54 5
Ruby
<lang ruby>require 'bigdecimal'
- square root of 2
sqrt2 = Object.new def sqrt2.a(n); n == 1 ? 1 : 2; end def sqrt2.b(n); 1; end
- Napier's constant
napier = Object.new def napier.a(n); n == 1 ? 2 : n - 1; end def napier.b(n); n == 1 ? 1 : n - 1; end
pi = Object.new def pi.a(n); n == 1 ? 3 : 6; end def pi.b(n); (2*n - 1)**2; end
- Estimates the value of a continued fraction _cfrac_, to _prec_
- decimal digits of precision. Returns a BigDecimal. _cfrac_ must
- respond to _cfrac.a(n)_ and _cfrac.b(n)_ for integer _n_ >= 1.
def estimate(cfrac, prec)
last_result = nil terms = prec
loop do # Estimate continued fraction for _n_ from 1 to _terms_. result = cfrac.a(terms) (terms - 1).downto(1) do |n| a = BigDecimal cfrac.a(n) b = BigDecimal cfrac.b(n) digits = [b.div(result, 1).exponent + prec, 1].max result = a + b.div(result, digits) end result = result.round(prec)
if result == last_result return result else # Double _terms_ and try again. last_result = result terms *= 2 end end
end
puts estimate(sqrt2, 50).to_s('F') puts estimate(napier, 50).to_s('F') puts estimate(pi, 10).to_s('F')</lang>
- Output:
$ ruby cfrac.rb 1.41421356237309504880168872420969807856967187537695 2.71828182845904523536028747135266249775724709369996 3.1415926536
Scala
Note that Scala-BigDecimal provides a precision of 34 digits. Therefore we take a limitation of 32 digits to avoiding rounding problems. <lang Scala>object CF extends App {
import Stream._ val sqrt2 = 1 #:: from(2,0) zip from(1,0) val napier = 2 #:: from(1) zip (1 #:: from(1)) val pi = 3 #:: from(6,0) zip (from(1,2) map {x=>x*x})
// reference values, source: wikipedia val refPi = "3.14159265358979323846264338327950288419716939937510" val refNapier = "2.71828182845904523536028747135266249775724709369995" val refSQRT2 = "1.41421356237309504880168872420969807856967187537694"
def calc(cf: Stream[(Int, Int)], numberOfIters: Int=200): BigDecimal = { (cf take numberOfIters toList).foldRight[BigDecimal](1)((a, z) => a._1+a._2/z) } def approx(cfV: BigDecimal, cfRefV: String): String = { val p: Pair[Char,Char] => Boolean = pair =>(pair._1==pair._2) ((cfV.toString+" "*34).substring(0,34) zip cfRefV.toString.substring(0,34)) .takeWhile(p).foldRight[String]("")((a:Pair[Char,Char],z)=>a._1+z) }
List(("sqrt2",sqrt2,50,refSQRT2),("napier",napier,50,refNapier),("pi",pi,3000,refPi)) foreach {t=> val (name,cf,iters,refV) = t val cfV = calc(cf,iters) println(name+":") println("ref value: "+refV.substring(0,34)) println("cf value: "+(cfV.toString+" "*34).substring(0,34)) println("precision: "+approx(cfV,refV)) println() }
}</lang> Output:
sqrt2: ref value: 1.41421356237309504880168872420969 cf value: 1.41421356237309504880168872420969 precision: 1.41421356237309504880168872420969 napier: ref value: 2.71828182845904523536028747135266 cf value: 2.71828182845904523536028747135266 precision: 2.71828182845904523536028747135266 pi: ref value: 3.14159265358979323846264338327950 cf value: 3.14159265358052780404906362935452 precision: 3.14159265358
For higher accuracy of pi we have to take more iterations. Unfortunately the foldRight function in calc isn't tail recursiv - therefore a stack overflow exception will be thrown for higher numbers of iteration, thus we have to implement an iterative way for calculation: <lang Scala>object CFI extends App {
import Stream._ val sqrt2 = 1 #:: from(2,0) zip from(1,0) val napier = 2 #:: from(1) zip (1 #:: from(1)) val pi = 3 #:: from(6,0) zip (from(1,2) map {x=>x*x})
// reference values, source: wikipedia val refPi = "3.14159265358979323846264338327950288419716939937510" val refNapier = "2.71828182845904523536028747135266249775724709369995" val refSQRT2 = "1.41421356237309504880168872420969807856967187537694"
def calc_i(cf: Stream[(Int, Int)], numberOfIters: Int=50): BigDecimal = { val cfl = cf take numberOfIters toList var z: BigDecimal = 1.0 for (i <- 0 to cfl.size-1 reverse) z=cfl(i)._1+cfl(i)._2/z z } def approx(cfV: BigDecimal, cfRefV: String): String = { val p: Pair[Char,Char] => Boolean = pair =>(pair._1==pair._2) ((cfV.toString+" "*34).substring(0,34) zip cfRefV.toString.substring(0,34)) .takeWhile(p).foldRight[String]("")((a:Pair[Char,Char],z)=>a._1+z) }
List(("sqrt2",sqrt2,50,refSQRT2),("napier",napier,50,refNapier),("pi",pi,50000,refPi)) foreach {t=> val (name,cf,iters,refV) = t val cfV = calc_i(cf,iters) println(name+":") println("ref value: "+refV.substring(0,34)) println("cf value: "+(cfV.toString+" "*34).substring(0,34)) println("precision: "+approx(cfV,refV)) println() }
}</lang> Output:
sqrt2: ref value: 1.41421356237309504880168872420969 cf value: 1.41421356237309504880168872420969 precision: 1.41421356237309504880168872420969 napier: ref value: 2.71828182845904523536028747135266 cf value: 2.71828182845904523536028747135266 precision: 2.71828182845904523536028747135266 pi: ref value: 3.14159265358979323846264338327950 cf value: 3.14159265358983426214354599901745 precision: 3.141592653589
Tcl
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
XPL0
The number of iterations (N) needed to get the 13 digits of accuracy was determined by experiment. <lang XPL0>include c:\cxpl\codes; int N; real A, B, F; [Format(1, 15); A:= 2.0; B:= 1.0; N:= 16; IntOut(0, N); CrLf(0); F:= 0.0; while N>=1 do [F:= B/(A+F); N:= N-1]; RlOut(0, 1.0+F); CrLf(0); RlOut(0, sqrt(2.0)); CrLf(0);
N:= 13; IntOut(0, N); CrLf(0); F:= 0.0; while N>=2 do [F:= float(N-1)/(float(N)+F); N:= N-1]; RlOut(0, 2.0 + 1.0/(1.0+F)); CrLf(0); RlOut(0, Exp(1.0)); CrLf(0);
N:= 10000; IntOut(0, N); CrLf(0); F:= 0.0; while N>=1 do [F:= float(sq(2*N-1))/(6.0+F); N:= N-1]; RlOut(0, 3.0+F); CrLf(0); RlOut(0, ACos(-1.0)); CrLf(0); ]</lang>
Output:
16 1.414213562372820 1.414213562373100 13 2.718281828459380 2.718281828459050 10000 3.141592653589540 3.141592653589790