Continued fraction
A number may be represented as a continued fraction (see Mathworld for more information) as follows:
You are encouraged to solve this task according to the task description, using any language you may know.
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
<lang haskell>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 terms cf2rat n = foldr (\(a,b) f -> (a%1) + ((b%1) / f)) (1%1) . take n
-- truncate after error is at most 1/p cf2rat_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</lang>
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>
NetRexx
<lang 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
- /
options replace format comments java crossref savelog symbols
Numeric Digits 30 Sqrt2 =1 napier=2 pi =3 a =0 b =0 Say 'SQRT2='.left(7) calc(sqrt2, 200,a,b,sqrt2,napier,pi) Say 'NAPIER='.left(7) calc(napier, 200,a,b,sqrt2,napier,pi) Say 'PI='.left(7) calc(pi, 200,a,b,sqrt2,napier,pi) Exit
method get_Coeffs(form,n,a,b,sqrt2,napier,pi) 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 a b
method calc(form,n,a,b,sqrt2,napier,pi) public static
temp=0 loop ni = n to 1 by -1 res=Get_Coeffs(form,ni,a,b,sqrt2,napier,pi) Parse res a b temp = b/(a + temp) end res=Get_Coeffs(form,0,a,b,sqrt2,napier,pi) Parse res a b return (A + temp)</lang>
Who could help me make a,b,sqrt2,napier,pi global (public) variables? This would simplify the solution:-)
PARI/GP
Partial solution for simple continued fractions. <lang parigp>back(v)=my(t=contfracpnqn(v));t[1,1]/t[2,1]*1. back(vector(100,i,2-(i==1)))</lang>
Output:
%1 = 1.4142135623730950488016887242096980786
PL/I
<lang 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; </lang> Result:
1.41421356237309505E+0000
Version for NAPIER: <lang>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;</lang>
2.71828182845904524E+0000
Version for SQRT2, NAPIER, PI <lang>/* 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; </lang> Output:
SQRT2= 1.41421356237309505 NAPIER= 2.71828182845904524 PI= 3.14159265358979349
Prolog
<lang 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 fractions
continued_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 2 sqrt_2_ab(0, 1, 1). sqrt_2_ab(_, 2, 1).
% definitions for napier napier_ab(0, 2, _). napier_ab(1, 1, 1). napier_ab(N, N, V) :- V is N - 1.
% definitions for pi pi_ab(0, 3, _). pi_ab(N, 6, V) :- V is (2 * N - 1)*(2 * N - 1). </lang> Output :
?- continued_fraction. sqrt(2) = 1.4142135623730951 e = 2.7182818284590455 pi = 3.141592622804847 true .
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
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.
<lang rexx>/*REXX program calculates & displays values of some continued fractions.*/
numeric digits 100 /*use 100 digits for the display.*/
terms=500 /*use 500 terms for calculations.*/
/*══════════════════════════════════════════════════════════════════════*/
a=1 rep(2); call tell '√2', cf(a)
/*══════════════════════════════════════════════════════════════════════*/
a=1 rep(1 2); call tell '√3', cf(a) /*also: 2∙sin(π/3) */
/*══════════════════════════════════════════════════════════════════════*/
/* ___ */
do N=2 to 17 /*generalized √ N */ a=1 rep(2); b=rep(N-1); call tell 'gen √'N, cf(a,b) end /*N*/
N=1/2; a=1 rep(2); b=rep(N-1); call tell 'gen √½', cf(a,b) /*══════════════════════════════════════════════════════════════════════*/
do j=1 for terms; a=a j; end; call tell 'e', cf(2 a, 1 a)
/*══════════════════════════════════════════════════════════════════════*/ a=copies(' 1',terms); call tell 'φ, phi', cf(1 a) /*══════════════════════════════════════════════════════════════════════*/
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] */
/*══════════════════════════════════════════════════════════════════════*/ terms=100000 /*use 100,000 terms for π calc.*/ a=rep(6)
do j=1 for terms; b=b (2*j-1)**2; end; call tell 'π, pi', cf(3 a,b)
exit /*stick a fork in it, we're done.*/ /*────────────────────────────────CF subroutine─────────────────────────*/ cf: procedure; parse arg C x,y; !=0; numeric digits digits()*2
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
end /*k*/
return !+C /*────────────────────────────────DIVZERO subroutine────────────────────*/ divZero: say; say '***error!***'; say 'division by zero.'; say; exit 13 /*────────────────────────────────REP subroutine────────────────────────*/ rep: parse arg rep; return space(copies(' 'rep,terms%words(rep))) /*────────────────────────────────TELL subroutine──────────────────────────────────────────*/ tell: parse arg ?,v; vv = left(format(v) / 1, 1 + digits())
say right(?,7) '=' vv ' α terms='left(a,50)
if symbol('B')=='VAR' then if b\== then say right(,7+2+digits()+2) ' ß terms='left(b,50) a=; b=; return</lang> output
√2 = 1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573 α 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
√3 = 1.732050807568877293527446341505872366942805253810380628055806979451933016908800037081146186757248576 α 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
gen √2 = 1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573 α 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
ß 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
gen √3 = 1.732050807568877293527446341505872366942805253810380628055806979451933016908800037081146186757248576 α 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
ß 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
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
ß 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
gen √5 = 2.236067977499789696409173668731276235440618359611525724270897245410520925637804899414414408378782275 α 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
ß 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
gen √6 = 2.449489742783178098197284074705891391965947480656670128432692567250960377457315026539859433104640235 α 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
ß 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
gen √7 = 2.645751311064590590501615753639260425710259183082450180368334459201068823230283627760392886474543611 α 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
ß 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
gen √8 = 2.828427124746190097603377448419396157139343750753896146353359475981464956924214077700775068655283145 α 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
ß 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
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
ß 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
gen √10 = 3.162277660168379331998893544432718533719555139325216826857504852792594438639238221344248108379300295 α 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
ß 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
gen √11 = 3.316624790355399849114932736670686683927088545589353597058682146116484642609043846708843399128290651 α 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
ß terms=10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
gen √12 = 3.464101615137754587054892683011744733885610507620761256111613958903866033817600074162292373514497151 α 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
ß terms=11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11
gen √13 = 3.605551275463989293119221267470495946251296573845246212710453056227166948293010445204619082018490718 α 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
ß terms=12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12
gen √14 = 3.741657386773941385583748732316549301756019807778726946303745467320035156306939027976809895194379572 α 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
ß terms=13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13 13
gen √15 = 3.872983346207416885179265399782399610832921705291590826587573766113483091936979033519287376858673518 α 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
ß terms=14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14
gen √16 = 4 α 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
ß terms=15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15 15
gen √17 = 4.123105625617660549821409855974077025147199225373620434398633573094954346337621593587863650810684297 α 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
ß terms=16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16
gen √½ = 0.707106781186547524400844362104849039284835937688474036588339868995366239231053519425193767163820786 α 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
ß terms=-0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5
e = 2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427 α terms= 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 2
φ, phi = 1.618033988749894848204586834365638117720309179805762862135448622705260462818902449707207204189391137 α 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
√3 = 1.732050807568877293527446341505872366942805253810380628055806979451933016908800037081146186757248576 α terms= 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
tan(1) = 1.557407724654902230506974807458360173087250772381520038383946605698861397151727289555099965202242984 α terms= 1 1 3 1 5 1 7 1 9 1 11 1 13 1 15 1 17 1 19 1 21 1
coth(1) = 1.313035285499331303636161246930847832912013941240452655543152967567084270461874382674679241480856303 α terms= 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37
coth(½) = 2.163953413738652848770004010218023117093738602150792272533574119296087634783339486574409418809750115 α terms= 6 10 14 18 22 26 30 34 38 42 46 50 54 58 62 66 70
π, pi = 3.141592653589792988470143264530440384041017830472772036746332303472711537960073664096818977224037083 α 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
ß terms= 1 9 25 49 81 121 169 225 289 361 441 529 625 729
Note: even with 200 digit accuracy and 100,000 terms, the calculate of π (pi) is only accurate to 15 digits.
Version 2 derived from PL/I
<lang rexx>/* 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)</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|
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 doubles.
<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