Continued fraction

From Rosetta Code
(Redirected from Continued fractions)
Jump to: navigation, search
Task
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 a0 = 1 then aN = 2. bN is always 1.

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

For Napier's Constant, use a0 = 2, then aN = N. b1 = 1 then bN = N − 1.

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

For Pi, use a0 = 3 then aN = 6. bN = (2N − 1)2.

\pi = 3 + \cfrac{1}{6 + \cfrac{9}{6 + \cfrac{25}{6 + \ddots}}}
See also:

Contents

[edit] Ada

(The source text for these examples can also be found on Bitbucket.)

Generic function for estimating continued fractions:

generic
type Scalar is digits <>;
 
with function A (N : in Natural) return Natural;
with function B (N : in Positive) return Natural;
function Continued_Fraction (Steps : in Natural) return Scalar;
function Continued_Fraction (Steps : in Natural) return Scalar is
function A (N : in Natural) return Scalar is (Scalar (Natural'(A (N))));
function B (N : in Positive) return Scalar is (Scalar (Natural'(B (N))));
 
Fraction : Scalar := 0.0;
begin
for N in reverse Natural range 1 .. Steps loop
Fraction := B (N) / (A (N) + Fraction);
end loop;
return A (0) + Fraction;
end Continued_Fraction;

Test program using the function above to estimate the square root of 2, Napiers constant and pi:

with Ada.Text_IO;
 
with Continued_Fraction;
 
procedure Test_Continued_Fractions is
type Scalar is digits 15;
 
package Square_Root_Of_2 is
function A (N : in Natural) return Natural is (if N = 0 then 1 else 2);
function B (N : in Positive) return Natural is (1);
 
function Estimate is new Continued_Fraction (Scalar, A, B);
end Square_Root_Of_2;
 
package Napiers_Constant is
function A (N : in Natural) return Natural is (if N = 0 then 2 else N);
function B (N : in Positive) return Natural is (if N = 1 then 1 else N-1);
 
function Estimate is new Continued_Fraction (Scalar, A, B);
end Napiers_Constant;
 
package Pi is
function A (N : in Natural) return Natural is (if N = 0 then 3 else 6);
function B (N : in Positive) return Natural is ((2 * N - 1) ** 2);
 
function Estimate is new Continued_Fraction (Scalar, A, B);
end Pi;
 
package Scalar_Text_IO is new Ada.Text_IO.Float_IO (Scalar);
use Ada.Text_IO, Scalar_Text_IO;
begin
Put (Square_Root_Of_2.Estimate (200), Exp => 0); New_Line;
Put (Napiers_Constant.Estimate (200), Exp => 0); New_Line;
Put (Pi.Estimate (10000), Exp => 0); New_Line;
end Test_Continued_Fractions;

[edit] Using only Ada 95 features

This example is exactly the same as the preceding one, but implemented using only Ada 95 features.

generic
type Scalar is digits <>;
 
with function A (N : in Natural) return Natural;
with function B (N : in Positive) return Natural;
function Continued_Fraction_Ada95 (Steps : in Natural) return Scalar;
function Continued_Fraction_Ada95 (Steps : in Natural) return Scalar is
function A (N : in Natural) return Scalar is
begin
return Scalar (Natural'(A (N)));
end A;
 
function B (N : in Positive) return Scalar is
begin
return Scalar (Natural'(B (N)));
end B;
 
Fraction : Scalar := 0.0;
begin
for N in reverse Natural range 1 .. Steps loop
Fraction := B (N) / (A (N) + Fraction);
end loop;
return A (0) + Fraction;
end Continued_Fraction_Ada95;
with Ada.Text_IO;
 
with Continued_Fraction_Ada95;
 
procedure Test_Continued_Fractions_Ada95 is
type Scalar is digits 15;
 
package Square_Root_Of_2 is
function A (N : in Natural) return Natural;
function B (N : in Positive) return Natural;
 
function Estimate is new Continued_Fraction_Ada95 (Scalar, A, B);
end Square_Root_Of_2;
 
package body Square_Root_Of_2 is
function A (N : in Natural) return Natural is
begin
if N = 0 then
return 1;
else
return 2;
end if;
end A;
 
function B (N : in Positive) return Natural is
begin
return 1;
end B;
end Square_Root_Of_2;
 
package Napiers_Constant is
function A (N : in Natural) return Natural;
function B (N : in Positive) return Natural;
 
function Estimate is new Continued_Fraction_Ada95 (Scalar, A, B);
end Napiers_Constant;
 
package body Napiers_Constant is
function A (N : in Natural) return Natural is
begin
if N = 0 then
return 2;
else
return N;
end if;
end A;
 
function B (N : in Positive) return Natural is
begin
if N = 1 then
return 1;
else
return N - 1;
end if;
end B;
end Napiers_Constant;
 
package Pi is
function A (N : in Natural) return Natural;
function B (N : in Positive) return Natural;
 
function Estimate is new Continued_Fraction_Ada95 (Scalar, A, B);
end Pi;
 
package body Pi is
function A (N : in Natural) return Natural is
begin
if N = 0 then
return 3;
else
return 6;
end if;
end A;
 
function B (N : in Positive) return Natural is
begin
return (2 * N - 1) ** 2;
end B;
end Pi;
 
package Scalar_Text_IO is new Ada.Text_IO.Float_IO (Scalar);
use Ada.Text_IO, Scalar_Text_IO;
begin
Put (Square_Root_Of_2.Estimate (200), Exp => 0); New_Line;
Put (Napiers_Constant.Estimate (200), Exp => 0); New_Line;
Put (Pi.Estimate (10000), Exp => 0); New_Line;
end Test_Continued_Fractions_Ada95;
Output:
 1.41421356237310
 2.71828182845905
 3.14159265358954

[edit] ATS

A fairly direct translation of the C version without using advanced features of the type system:

(* a coefficient function creates double values from in paramters *)
typedef coeff_f = int -> double
 
(* a continued fraction is described by a record of two coefficent functions a and b *)
typedef frac = @{a= coeff_f, b= coeff_f}
 
(* recursive definition of the approximation *)
fun calc_rec(f: frac, n: int, r: double): double =
if n = 0 then
f.a(0) + r // base case
else let // recursive case
val a: double = f.a(n) // record access
val b: double = f.b(n)
val s: double = b / (a + r)
//val () = printf("r%d = %f, a = %f, b = %f, r%d = %f\n", @(n, r, a, b, (n-1), s))
in calc_rec(f, n - 1, s) end
 
fun calc(f: frac, n: int): double =
calc_rec(f, n, 0.0)
 
(* anonymous coefficient functions created with lam (short for lambda) *)
val sqrt2 = @{ a= lam (n: int): double => if n = 0 then 1.0 else 2.0,
b= lam (n: int): double => 1.0 }
 
val napier = @{ a= lam (n: int): double => if n = 0 then 2.0 else 1.0 * (n),
b= lam (n: int): double => if n = 1 then 1.0 else n - 1.0 }
 
fun square(x: double): double = x * x
 
val pi = @{ a= lam (n: int): double => if n = 0 then 3.0 else 6.0,
b= lam (n: int): double => square (2.0 * n - 1) }
 
implement main () = begin
printf ("sqrt2 = %f\n", @(calc(sqrt2, 100)));
printf ("napier = %f\n", @(calc(napier, 100)));
printf (" pi = %f\n", @(calc( pi , 100)));
end

[edit] Axiom

Axiom provides a ContinuedFraction domain:

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
Output:
   (1)  1.4142135623 730950488
Type: Float
 
(2) 2.7182818284 590452354
Type: Float
 
(3) 3.1415926538 39792926
Type: Float

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:

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

[edit] BBC BASIC

      *FLOAT64
@% = &1001010
 
PRINT "SQR(2) = " ; FNcontfrac(1, 1, "2", "1")
PRINT " e = " ; FNcontfrac(2, 1, "N", "N")
PRINT " PI = " ; FNcontfrac(3, 1, "6", "(2*N+1)^2")
END
 
REM a$ and b$ are functions of N
DEF FNcontfrac(a0, b1, a$, b$)
LOCAL N, expr$
REPEAT
N += 1
expr$ += STR$(EVAL(a$)) + "+" + STR$(EVAL(b$)) + "/("
UNTIL LEN(expr$) > (65500 - N)
= a0 + b1 / EVAL (expr$ + "1" + STRING$(N, ")"))
Output:
SQR(2) = 1.414213562373095
     e = 2.718281828459046
    PI = 3.141592653588017

[edit] C

Works with: ANSI C
/* calculate approximations for continued fractions */
#include <stdio.h>
 
/* kind of function that returns a series of coefficients */
typedef double (*coeff_func)(unsigned n);
 
/* calculates the specified number of expansions of the continued fraction
* described by the coefficient series f_a and f_b */

double calc(coeff_func f_a, coeff_func f_b, unsigned expansions)
{
double a, b, r;
a = b = r = 0.0;
 
unsigned i;
for (i = expansions; i > 0; i--) {
a = f_a(i);
b = f_b(i);
r = b / (a + r);
}
a = f_a(0);
 
return a + r;
}
 
/* series for sqrt(2) */
double sqrt2_a(unsigned n)
{
return n ? 2.0 : 1.0;
}
 
double sqrt2_b(unsigned n)
{
return 1.0;
}
 
/* series for the napier constant */
double napier_a(unsigned n)
{
return n ? n : 2.0;
}
 
double napier_b(unsigned n)
{
return n > 1.0 ? n - 1.0 : 1.0;
}
 
/* series for pi */
double pi_a(unsigned n)
{
return n ? 6.0 : 3.0;
}
 
double pi_b(unsigned n)
{
double c = 2.0 * n - 1.0;
 
return c * c;
}
 
int main(void)
{
double sqrt2, napier, pi;
 
sqrt2 = calc(sqrt2_a, sqrt2_b, 1000);
napier = calc(napier_a, napier_b, 1000);
pi = calc(pi_a, pi_b, 1000);
 
printf("%12.10g\n%12.10g\n%12.10g\n", sqrt2, napier, pi);
 
return 0;
}
Output:
 1.414213562
 2.718281828
 3.141592653

[edit] C++

#include <iomanip>
#include <iostream>
#include <tuple>
 
typedef std::tuple<double,double> coeff_t; // coefficients type
typedef coeff_t (*func_t)(int); // callback function type
 
double calc(func_t func, int n)
{
double a, b, temp = 0;
for (; n > 0; --n) {
std::tie(a, b) = func(n);
temp = b / (a + temp);
}
std::tie(a, b) = func(0);
return a + temp;
}
 
coeff_t sqrt2(int n)
{
return coeff_t(n > 0 ? 2 : 1, 1);
}
 
coeff_t napier(int n)
{
return coeff_t(n > 0 ? n : 2, n > 1 ? n - 1 : 1);
}
 
coeff_t pi(int n)
{
return coeff_t(n > 0 ? 6 : 3, (2 * n - 1) * (2 * n - 1));
}
 
int main()
{
std::streamsize old_prec = std::cout.precision(15); // set output digits
std::cout
<< calc(sqrt2, 20) << '\n'
<< calc(napier, 15) << '\n'
<< calc(pi, 10000) << '\n'
<< std::setprecision(old_prec); // reset precision
}
Output:
1.41421356237309
2.71828182845905
3.14159265358954

[edit] 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)
Output:
> coffee continued_fraction.coffee 
1.4142135623730951
1.4142135623730951
---
2.718281828459045
2.7182818284590455
---
3.141592653589793
3.141592653589793
---

[edit] Common Lisp

Translation of: C++
(defun estimate-continued-fraction (generator n)
(let ((temp 0))
(loop for n1 from n downto 1
do (multiple-value-bind (a b)
(funcall generator n1)
(setf temp (/ b (+ a temp)))))
(+ (funcall generator 0) temp)))
 
(format t "sqrt(2) = ~a~%" (coerce (estimate-continued-fraction
(lambda (n)
(values (if (> n 0) 2 1) 1)) 20)
'double-float))
(format t "napier's = ~a~%" (coerce (estimate-continued-fraction
(lambda (n)
(values (if (> n 0) n 2)
(if (> n 1) (1- n) 1))) 15)
'double-float))
 
(format t "pi = ~a~%" (coerce (estimate-continued-fraction
(lambda (n)
(values (if (> n 0) 6 3)
(* (1- (* 2 n))
(1- (* 2 n))))) 10000)
'double-float))
Output:
sqrt(2) = 1.4142135623730947d0
napier's = 2.7182818284590464d0
pi = 3.141592653589543d0

[edit] Chapel

Functions don't take other functions as arguments, so I wrapped them in a dummy record each.

proc calc(f, n) {
var r = 0.0;
 
for k in 1..n by -1 {
var v = f.pair(k);
r = v(2) / (v(1) + r);
}
 
return f.pair(0)(1) + r;
}
 
record Sqrt2 {
proc pair(n) {
return (if n == 0 then 1 else 2,
1);
}
}
 
record Napier {
proc pair(n) {
return (if n == 0 then 2 else n,
if n == 1 then 1 else n - 1);
}
}
record Pi {
proc pair(n) {
return (if n == 0 then 3 else 6,
(2*n - 1)**2);
}
}
 
config const n = 200;
writeln(calc(new Sqrt2(), n));
writeln(calc(new Napier(), n));
writeln(calc(new Pi(), n));

[edit] D

import std.stdio, std.functional, std.traits;
 
FP calc(FP, F)(in F fun, in int n) pure nothrow if (isCallable!F) {
FP temp = 0;
 
foreach_reverse (immutable ni; 1 .. n + 1) {
immutable p = fun(ni);
temp = p[1] / (FP(p[0]) + temp);
}
return fun(0)[0] + temp;
}
 
int[2] fSqrt2(in int n) pure nothrow {
return [n > 0 ? 2 : 1, 1];
}
 
int[2] fNapier(in int n) pure nothrow {
return [n > 0 ? n : 2, n > 1 ? (n - 1) : 1];
}
 
int[2] fPi(in int n) pure nothrow {
return [n > 0 ? 6 : 3, (2 * n - 1) ^^ 2];
}
 
alias print = curry!(writefln, "%.19f");
 
void main() {
calc!real(&fSqrt2, 200).print;
calc!real(&fNapier, 200).print;
calc!real(&fPi, 200).print;
}
Output:
1.4142135623730950487
2.7182818284590452354
3.1415926228048469486

[edit] Erlang

 
-module(continued).
-compile([export_all]).
 
pi_a (0) -> 3;
pi_a (_N) -> 6.
 
pi_b (N) ->
(2*N-1)*(2*N-1).
 
sqrt2_a (0) ->
1;
sqrt2_a (_N) ->
2.
 
sqrt2_b (_N) ->
1.
 
nappier_a (0) ->
2;
nappier_a (N) ->
N.
 
nappier_b (1) ->
1;
nappier_b (N) ->
N-1.
 
continued_fraction(FA,_FB,0) -> FA(0);
continued_fraction(FA,FB,N) ->
continued_fraction(FA,FB,N-1,FB(N)/FA(N)).
 
continued_fraction(FA,_FB,0,Acc) -> FA(0) + Acc;
continued_fraction(FA,FB,N,Acc) ->
continued_fraction(FA,FB,N-1,FB(N)/ (FA(N) + Acc)).
 
test_pi (N) ->
continued_fraction(fun pi_a/1,fun pi_b/1,N).
 
test_sqrt2 (N) ->
continued_fraction(fun sqrt2_a/1,fun sqrt2_b/1,N).
 
test_nappier (N) ->
continued_fraction(fun nappier_a/1,fun nappier_b/1,N).
 
Output:
 
29> continued:test_pi(1000).
3.141592653340542
30> continued:test_sqrt2(1000).
1.4142135623730951
31> continued:test_nappier(1000).
2.7182818284590455
 

[edit] 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.

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
Output:
 Square root of 2: 1.414213562373095048801688724209
Napier's constant: 2.718281828459045235360287471352
               Pi: 3.1415926538

[edit] Forth

Translation of: D
: 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
;
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

[edit] Fortran

module continued_fractions
implicit none
 
integer, parameter :: long = selected_real_kind(7,99)
 
type continued_fraction
integer :: a0, b1
procedure(series), pointer, nopass :: a, b
end type
 
interface
pure function series (n)
integer, intent(in) :: n
integer :: series
end function
end interface
 
contains
 
pure function define_cf (a0,a,b1,b) result(x)
integer, intent(in) :: a0
procedure(series) :: a
integer, intent(in), optional :: b1
procedure(series), optional :: b
type(continued_fraction) :: x
x%a0 = a0
x%a => a
if ( present(b1) ) then
x%b1 = b1
else
x%b1 = 1
end if
if ( present(b) ) then
x%b => b
else
x%b => const_1
end if
end function define_cf
 
pure integer function const_1(n)
integer,intent(in) :: n
const_1 = 1
end function
 
pure real(kind=long) function output(x,iterations)
type(continued_fraction), intent(in) :: x
integer, intent(in) :: iterations
integer :: i
output = x%a(iterations)
do i = iterations-1,1,-1
output = x%a(i) + (x%b(i+1) / output)
end do
output = x%a0 + (x%b1 / output)
end function output
 
end module continued_fractions
 
 
program examples
use continued_fractions
 
type(continued_fraction) :: sqr2,napier,pi
 
sqr2 = define_cf(1,a_sqr2)
napier = define_cf(2,a_napier,1,b_napier)
pi = define_cf(3,a=a_pi,b=b_pi)
 
write (*,*) output(sqr2,10000)
write (*,*) output(napier,10000)
write (*,*) output(pi,10000)
 
contains
 
pure integer function a_sqr2(n)
integer,intent(in) :: n
a_sqr2 = 2
end function
 
pure integer function a_napier(n)
integer,intent(in) :: n
a_napier = n
end function
 
pure integer function b_napier(n)
integer,intent(in) :: n
b_napier = n-1
end function
 
pure integer function a_pi(n)
integer,intent(in) :: n
a_pi = 6
end function
 
pure integer function b_pi(n)
integer,intent(in) :: n
b_pi = (2*n-1)*(2*n-1)
end function
 
end program examples
Output:
   1.4142135623730951
   2.7182818284590455
   3.1415926535895435

[edit] 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())
}
Output:
sqrt2: 1.4142135623730965
nap:   2.7182818284590455
pi:    3.141623806667839

[edit] 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]
Output:
1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572735013846230912297024924836055850737212644121497099935831413222665927505592755799950501152782060571
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019
3.141592653297590947683406834261190738869139611505752231394089152890909495973464508817163306557131591579057202097715021166512662872910519439747609829479577279606075707015622200744006783543589980682386
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

[edit] 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

[edit] Mathematica

sqrt2=Function[n,{1,Transpose@{Array[2&,n],Array[1&,n]}}];
napier=Function[n,{2,Transpose@{Range[n],Prepend[Range[n-1],1]}}];
pi=Function[n,{3,Transpose@{Array[6&,n],Array[(2#-1)^2&,n]}}];
approx=Function[l,
N[Divide@@First@Fold[{{#2.#[[;;,1]],#2.#[[;;,2]]},#[[1]]}&,{{l[[2,1,1]]l[[1]]+l[[2,1,2]],l[[2,1,1]]},{l[[1]],1}},l[[2,2;;]]],10]];
r2=approx/@{sqrt2@#,napier@#,pi@#}&@10000;r2//TableForm
Output:
1.414213562
2.718281828
3.141592654

[edit] 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 */

[edit] 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
* 08.09.2012 Walter Pachl simplified (with the help of a friend)
**********************************************************************/

options replace format comments java crossref savelog symbols
class CFB public
 
properties static
Numeric Digits 30
Sqrt2 =1
napier=2
pi =3
a =0
b =0
 
method main(args = String[]) public static
Say 'SQRT2='.left(7) calc(sqrt2, 200)
Say 'NAPIER='.left(7) calc(napier, 200)
Say 'PI='.left(7) calc(pi, 200)
Return
 
method get_Coeffs(form,n) 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
 
method calc(form,n) public static
temp=0
loop ni = n to 1 by -1
Get_Coeffs(form,ni)
temp = b/(a + temp)
end
Get_Coeffs(form,0)
return (a + temp)

Who could help me make a,b,sqrt2,napier,pi global (public) variables? This would simplify the solution:-)

I got this help and simplified the program.

However, I am told that 'my' value of pi is incorrect. I will investigate!

Apparently the coefficients given in the task description are only good for an approximation. One should, therefore, not SHOW more that 15 digits. See http://de.wikipedia.org/wiki/Kreiszahl

See Rexx for a better computation

[edit] OCaml

let pi = 3, fun n -> ((2*n-1)*(2*n-1), 6)
and nap = 2, fun n -> (max 1 (n-1), n)
and root2 = 1, fun n -> (1, 2) in
 
let eval (i,f) k =
let rec frac n =
let a, b = f n in
float a /. (float b +.
if n >= k then 0.0 else frac (n+1)) in
float i +. frac 1 in
 
Printf.printf "sqrt(2)\t= %.15f\n" (eval root2 1000);
Printf.printf "e\t= %.15f\n" (eval nap 1000);
Printf.printf "pi\t= %.15f\n" (eval pi 1000);

Output (inaccurate due to too few terms):

sqrt(2)	= 1.414213562373095
e	= 2.718281828459046
pi	= 3.141592653340542

[edit] PARI/GP

Partial solution for simple continued fractions.

back(v)=my(t=contfracpnqn(v));t[1,1]/t[2,1]*1.
back(vector(100,i,2-(i==1)))

Output:

%1 = 1.4142135623730950488016887242096980786

[edit] Perl

We'll use closures to implement the infinite lists of coeffficients.

sub continued_fraction {
my ($a, $b, $n) = (@_[0,1], $_[2] // 100);
 
$a->() + ($n && $b->() / continued_fraction($a, $b, $n-1));
}
 
printf "√2 ≈ %.9f\n", continued_fraction do { my $n; sub { $n++ ? 2 : 1 } }, sub { 1 };
printf "e ≈ %.9f\n", continued_fraction do { my $n; sub { $n++ || 2 } }, do { my $n; sub { $n++ || 1 } };
printf "π ≈ %.9f\n", continued_fraction do { my $n; sub { $n++ ? 6 : 3 } }, do { my $n; sub { (2*$n++ + 1)**2 } }, 1_000;
printf "π/2 ≈ %.9f\n", continued_fraction do { my $n; sub { 1/($n++ || 1) } }, sub { 1 }, 1_000;
Output:
√2  ≈ 1.414213562
e   ≈ 2.718281828
π   ≈ 3.141592653
π/2 ≈ 1.570717797

[edit] Perl 6

sub continued-fraction(:@a, :@b, Int :$n = 100)
{
my $x = @a[$n - 1];
$x = @a[$_ - 1] + @b[$_] / $x for reverse 1 ..^ $n;
$x;
}
 
printf "√2 ≈ %.9f\n", continued-fraction(:a(1, 2 xx *), :b(*, 1 xx *));
printf "e ≈ %.9f\n", continued-fraction(:a(2, 1 .. *), :b(*, 1, 1 .. *));
printf "π ≈ %.9f\n", continued-fraction(:a(3, 6 xx *), :b(*, [\+] 1, (8, 16 ... *)), :n(1000));
Output:
√2 ≈ 1.414213562
e  ≈ 2.718281828
π  ≈ 3.141592654

A more original and a bit more abstract method would consist in viewing a continued fraction on rank n as a function of a variable x:

\mathrm{CF}_3(x) = a_0 + \cfrac{b_1}{a_1 + \cfrac{b_2}{a_2 + \cfrac{b_3}{a_3 + x}}}

Or, more consistently:

\mathrm{CF}_3(x) = a_0 + \cfrac{b_0}{a_1 + \cfrac{b_1}{a_2 + \cfrac{b_2}{a_3 + \cfrac{b_3}{x}}}}

Viewed as such, CFn(x) could be written recursively:

\mathrm{CF}_n(x) = \mathrm{CF}_{n-1}(a_n + \frac{b_n}{x})

Or in other words:

\mathrm{CF}_n= \mathrm{CF}_{n-1}\circ f_n = \mathrm{CF}_{n-2}\circ f_{n-1}\circ f_n=\ldots=f_0\circ f_1 \ldots \circ f_n

where f_n(x) = a_n + \frac{b_n}{x}

Perl6 allows us to define a custom composition operator. We can then use it with the triangular reduction metaoperator, and map each resulting function with an infinite value for x (any value would do actually, but infinite make it consistent with this particular task).

sub infix:<>(&f, &g) { -> $x { &f(&g($x)) } }
sub continued-fraction(@a, @b, $x = Inf) {
map { .($x) },
[\⚬] map { @a[$_] + @b[$_] / * }, 0 .. *
}
 
printf "√2 ≈ %.9f\n", continued-fraction((1, 2 xx *), (1 xx *))[10];
printf "e ≈ %.9f\n", continued-fraction((2, 1 .. *), (1, 1 .. *))[10];
printf "π ≈ %.9f\n", continued-fraction((3, 6 xx *), ((1, 3, 5 ... *) X** 2))[100];

The main advantage is that the definition of the function does not need to know for which rank n it is computed. This is arguably closer to the mathematical definition.

[edit] 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;
Output:
 1.41421356237309505E+0000 

Version for NAPIER:

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;
 2.71828182845904524E+0000 

Version for SQRT2, NAPIER, PI

/* 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;
Output:
SQRT2=     1.41421356237309505
NAPIER=    2.71828182845904524
PI=        3.14159265358979349

[edit] 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).
Output:
 ?- continued_fraction.
sqrt(2) = 1.4142135623730951
e       = 2.7182818284590455
pi      = 3.141592622804847
true .

[edit] Python

Works with: Python version 2.6+ and 3.x
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

[edit] Fast iterative version

Translation of: D
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)
Output:
1.4142135623730950488016887242096980785696718753770
2.7182818284590452353602874713526624977572470937000
3.1415926839198062649342019294083175420335002640134

[edit] Racket

[edit] Using Doubles

This version uses standard double precision floating point numbers:

 
#lang racket
(define (calc cf n)
(match/values (cf 0)
[(a0 b0)
(+ a0
(for/fold ([t 0.0]) ([i (in-range (+ n 1) 0 -1)])
(match/values (cf i)
[(a b) (/ b (+ a t))])))]))
 
(define (cf-sqrt i) (values (if (> i 0) 2 1) 1))
(define (cf-napier i) (values (if (> i 0) i 2) (if (> i 1) (- i 1) 1)))
(define (cf-pi i) (values (if (> i 0) 6 3) (sqr (- (* 2 i) 1))))
 
(calc cf-sqrt 200)
(calc cf-napier 200)
(calc cf-pi 200)
 

Output:

 
1.4142135623730951
2.7182818284590455
3.1415926839198063
 

[edit] Version - Using Doubles

This versions uses big floats (arbitrary precision floating point):

 
#lang racket
(require math)
(bf-precision 2048) ; in bits
 
(define (calc cf n)
(match/values (cf 0)
[(a0 b0)
(bf+ (bf a0)
(for/fold ([t (bf 0)]) ([i (in-range (+ n 1) 0 -1)])
(match/values (cf i)
[(a b) (bf/ (bf b) (bf+ (bf a) t))])))]))
 
(define (cf-sqrt i) (values (if (> i 0) 2 1) 1))
(define (cf-napier i) (values (if (> i 0) i 2) (if (> i 1) (- i 1) 1)))
(define (cf-pi i) (values (if (> i 0) 6 3) (sqr (- (* 2 i) 1))))
 
(calc cf-sqrt 200)
(calc cf-napier 200)
(calc cf-pi 200)
 

Output:

 
(bf #e1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462309122970249248360558507372126441214970999358960036439214262599769155193770031712304888324413327207659690547583107739957489062466508437105234564161085482146113860092820802430986649987683947729823677905101453725898480737256099166805538057375451207262441039818826744940289448489312217214883459060818483750848688583833366310472320771259749181255428309841375829513581694269249380272698662595131575038315461736928338289219865139248048189188905788104310928762952913687232022557677738108337499350045588767581063729)
(bf #e2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193200305992181741359662904357290033429526059563073813232862794349076323382988075319525101901157383418793070215408914993488416750924476146066808226480016847741185374234544243710753907774499206955170276183860626133138458300075204493382656029760673711320070932870912744374704723624212700454495421842219077173525899689811474120614457405772696521446961165559468253835854362096088934714907384964847142748311021268578658461064714894910680584249490719358138073078291397044213736982988247857479512745588762993966446075)
(bf #e3.14159268391980626493420192940831754203350026401337226640663040854412059241988978103217808449508253393479795573626200366332733859609651462659489470805432281782785922056335606047700127154963266242144951481397480765182268219697420028007903565511884267297358842935537138583640066772149177226656227031792115896439889412205871076985598822285367358003457939603015797225018209619662200081521930463480571130673429337524564941105654923909951299948539893933654293161126559643573974163405197696633200469475250152247413175932572922175467223988860975105100904322239324381097207835036465269418118204894206705789759765527734394105147)
 

[edit] REXX

[edit] 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 (and also ½). In addition, ½π was calculated (as described in the talk page under Gold Credit).

More code is used for nicely formatting the output than the continued fraction calculation.

/*REXX program calculates and displays values of some specific continued*/
/*───────────── fractions (along with their α and ß terms). */
/*───────────── Continued fractions: also known as anthyphairetic ratio.*/
T=500 /*use 500 terms for calculations.*/
showDig=100; numeric digits 2*showDig /*use 100 digits for the display.*/
a=; @=; b= /*omitted ß terms are assumed = 1*/
/*══════════════════════════════════════════════════════════════════════*/
a=1 rep(2); call tell '√2'
/*══════════════════════════════════════════════════════════════════════*/
a=1 rep(1 2); call tell '√3' /*also: 2∙sin(π/3) */
/*══════════════════════════════════════ ___ ════════════════════════*/
/*generalized √ N */
do N=2 to 11; a=1 rep(2); b=rep(N-1); call tell 'gen √'N; end
N=1/2; a=1 rep(2); b=rep(N-1); call tell 'gen √½'
/*══════════════════════════════════════════════════════════════════════*/
do j=1 for T; a=a j; end; b=1 a; a=2 a; call tell 'e'
/*══════════════════════════════════════════════════════════════════════*/
do j=1 for T by 2; a=a j; b=b j+1; end; call tell '1÷[√e-1]'
/*══════════════════════════════════════════════════════════════════════*/
do j=1 for T; a=a j; end; b=a; a=0 a; call tell '1÷[e-1]'
/*══════════════════════════════════════════════════════════════════════*/
a=1 rep(1); call tell 'φ, phi'
/*══════════════════════════════════════════════════════════════════════*/
a=1; do j=1 for T by 2; a=a j 1; end; call tell 'tan(1)'
/*══════════════════════════════════════════════════════════════════════*/
a=1; do j=1 for T; a=a 2*j+1; end; call tell 'coth(1)'
/*══════════════════════════════════════════════════════════════════════*/
a=2; do j=1 for T; a=a 4*j+2; end; call tell 'coth(½)' /*also: [e+1] ÷ [e-1] */
/*══════════════════════════════════════════════════════════════════════*/
T=10000
a=1 rep(2)
do j=1 for T by 2; b=b j**2; end; call tell '4÷π'
/*══════════════════════════════════════════════════════════════════════*/
T=10000
a=1; do j=1 for T; a=a 1/j; @=@ '1/'j; end; call tell '½π, ½pi'
/*══════════════════════════════════════════════════════════════════════*/
T=10000
a=0 1 rep(2)
do j=1 for T by 2; b=b j**2; end; b=4 b; call tell 'π, pi'
/*══════════════════════════════════════════════════════════════════════*/
T=10000
a=0; do j=1 for T; a=a j*2-1; b=b j**2; end; b=4 b; call tell 'π, pi'
/*══════════════════════════════════════════════════════════════════════*/
T=100000
a=3 rep(6)
do j=1 for T by 2; b=b j**2; end; call tell 'π, pi'
exit /*stick a fork in it, we're done.*/
 
/*────────────────────────────────CF subroutine─────────────────────────*/
cf: procedure; parse arg C x,y;  !=0; numeric digits digits()+5
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 /*here's a binary mosh pit.*/
end /*k*/
return !+C
/*────────────────────────────────DIVZERO subroutine────────────────────*/
divZero: say; say '***error!***'; say 'division by zero.'; say; exit 13
/*────────────────────────────────GETT subroutine───────────────────────*/
getT: parse arg stuff,width,ma,mb,_
do m=1; mm=m+ma; mn=max(1,m-mb); w=word(stuff,m)
w=right(w,max(length(word(As,mm)),length(word(Bs,mn)),length(w)))
if length(_ w)>width then leave /*stop getting terms?*/
_=_ w /*whole, don't chop. */
end /*m*/ /*done building terms*/
return strip(_) /*strip leading blank*/
/*────────────────────────────────REP subroutine────────────────────────*/
rep: parse arg rep; return space(copies(' 'rep, T%words(rep)))
/*────────────────────────────────RF subroutine─────────────────────────*/
rf: parse arg xxx,z
do m=1 for T; w=word(xxx,m)  ; if w=='1/1' | w=1 then w=1
if w=='1/2' | w=1/2 then w='½'; if w=-.5 then w='-½'
if w=='1/4' | w=1/4 then w='¼'; if w=-.25 then w='-¼'
z=z w
end
return z /*done re-formatting.*/
/*────────────────────────────────TELL subroutine───────────────────────*/
tell: parse arg ?; v=cf(a,b); numeric digits showdig; As=rf(@ a); Bs=rf(b)
say right(?,8) '=' left(v/1,showdig) ' α terms= ' getT(As,72 ,0,1)
if b\=='' then say right('',8+2+showdig+1) ' ß terms= ' getT(Bs,72-2,1,0)
a=; @=; b=; return

output

      √2 = 1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157   α 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 2 2 2 2 2 2 2 2 2 2 2
      √3 = 1.73205080756887729352744634150587236694280525381038062805580697945193301690880003708114618675724857   α 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 1 2 1 2 1 2 1 2 1 2 1
  gen √2 = 1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157   α 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 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 1 1 1 1 1 1 1 1 1 1
  gen √3 = 1.73205080756887729352744634150587236694280525381038062805580697945193301690880003708114618675724857   α 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 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 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 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 3 3 3 3 3 3 3 3 3 3
  gen √5 = 2.23606797749978969640917366873127623544061835961152572427089724541052092563780489941441440837878227   α 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 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 4 4 4 4 4 4 4 4 4 4
  gen √6 = 2.44948974278317809819728407470589139196594748065667012843269256725096037745731502653985943310464023   α 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 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 5 5 5 5 5 5 5 5 5 5
  gen √7 = 2.64575131106459059050161575363926042571025918308245018036833445920106882323028362776039288647454361   α 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 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 6 6 6 6 6 6 6 6 6 6
  gen √8 = 2.82842712474619009760337744841939615713934375075389614635335947598146495692421407770077506865528314   α 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 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 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 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 8 8 8 8 8 8 8 8 8 8
 gen √10 = 3.16227766016837933199889354443271853371955513932521682685750485279259443863923822134424810837930029   α 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 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 9 9 9 9 9 9 9 9 9 9
 gen √11 = 3.31662479035539984911493273667068668392708854558935359705868214611648464260904384670884339912829065   α 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
                                                                                                                  ß terms=    10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
  gen √½ = 0.70710678118654752440084436210484903928483593768847403658833986899536623923105351942519376716382078   α 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
                                                                                                                  ß terms=    -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½ -½
       e = 2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642   α terms=  2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
                                                                                                                  ß terms=    1 1 2 3 4 5 6 7 8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
1÷[√e-1] = 1.54149408253679828413110344447251463834045923684188210947413695663754263914331480707182572408500774   α terms=  1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49
                                                                                                                  ß terms=    2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48
 1÷[e-1] = 0.58197670686932642438500200510901155854686930107539613626678705964804381739166974328720470940487505   α terms=  0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
                                                                                                                  ß terms=    1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
  φ, phi = 1.61803398874989484820458683436563811772030917980576286213544862270526046281890244970720720418939113   α 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 1 1 1 1 1 1 1 1 1 1 1
  tan(1) = 1.55740772465490223050697480745836017308725077238152003838394660569886139715172728955509996520224298   α terms=  1 1 1 3 1 5 1 7 1 9 1 11 1 13 1 15 1 17 1 19 1 21 1 23 1 25 1 27 1 29 1
 coth(1) = 1.31303528549933130363616124693084783291201394124045265554315296756708427046187438267467924148085630   α terms=  1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49
 coth(½) = 2.16395341373865284877000401021802311709373860215079227253357411929608763478333948657440941880975011   α terms=  2 6 10 14 18 22 26 30 34 38 42 46 50 54 58 62 66 70 74 78 82 86 90 94
     4÷π = 1.27283479368898552763028955877501991659132774562422849516943825241995907438793488819711543252100078   α terms=  1 2 2  2  2  2   2   2   2   2   2   2   2   2   2   2   2    2    2
                                                                                                                  ß terms=    1 9 25 49 81 121 169 225 289 361 441 529 625 729 841 961 1089 1225
 ½π, ½pi = 1.57071779679495409624134340842172803914665374867756685353559415360226497865248110814032818773478787   α terms=  1 ½ 1/3 ¼ 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17
   π, pi = 3.14259165433954305090112773725220456615353825631695587367530386050342717161955770321769607013860474   α terms=  0 1 2 2  2  2  2   2   2   2   2   2   2   2   2   2   2   2    2    2
                                                                                                                  ß terms=    4 1 9 25 49 81 121 169 225 289 361 441 529 625 729 841 961 1089 1225
   π, pi = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706   α terms=  0 1 3 5 7  9 11 13 15 17 19  21  23  25  27  29  31  33  35  37  39  41
                                                                                                                  ß terms=    4 1 4 9 16 25 36 49 64 81 100 121 144 169 196 225 256 289 324 361 400
   π, pi = 3.14159265358979298847014326453044038404101783047277203674633230347271153796007366409681897722403708   α terms=  3 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 841 961 1089 1225

Note: even with 200 digit accuracy and 100,000 terms, the last calculation of π (pi) is only accurate to 15 digits.

[edit] Version 2 derived from PL/I

/* 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)

[edit] Version 3 better approximation

/* REXX *************************************************************
* The task description specifies a continued fraction for pi
* that gives a reasonable approximation.
* Literature shows a better CF that yields pi with a precision of
* 200 digits.
* http://de.wikipedia.org/wiki/Kreiszahl
* 1
* pi = 3 + ------------------------
* 1
* 7 + --------------------
* 1
* 15 + ---------------
* 1
* 1 + -----------
*
* 292 + ...
*
* This program uses that CF and shows the first 50 digits
* PI =3.1415926535897932384626433832795028841971693993751...
* PIX=3.1415926535897932384626433832795028841971693993751...
* 201 correct digits
* 18.09.2012 Walter Pachl
**********************************************************************/

pi='3.1415926535897932384626433832795028841971'||,
'693993751058209749445923078164062862089986280348'||,
'253421170679821480865132823066470938446095505822'||,
'317253594081284811174502841027019385211055596446'||,
'229489549303819644288109756659334461284756482337'||,
'867831652712019091456485669234603486104543266482'||,
'133936072602491412737245870066063155881748815209'||,
'209628292540917153643678925903600113305305488204'||,
'665213841469519415116094330572703657595919530921'||,
'861173819326117931051185480744623799627495673518'||,
'857527248912279381830119491298336733624'
Numeric Digits 1000
al='7 15 1 292 1 1 1 2 1 3 1 14 2 1 1 2 2 2 2 1 84 2',
'1 1 15 3 13 1 4 2 6 6 99 1 2 2 6 3 5 1 1 6 8 1 7 1 2',
'3 7 1 2 1 1 12 1 1 1 3 1 1 8 1 1 2 1 6 1 1 5 2 2 3 1',
'2 4 4 16 1 161 45 1 22 1 2 2 1 4 1 2 24 1 2 1 3 1 2',
'1 1 10 2 5 4 1 2 2 8 1 5 2 2 26 1 4 1 1 8 2 42 2 1 7',
'3 3 1 1 7 2 4 9 7 2 3 1 57 1 18 1 9 19 1 2 18 1 3 7',
'30 1 1 1 3 3 3 1 2 8 1 1 2 1 15 1 2 13 1 2 1 4 1 12',
'1 1 3 3 28 1 10 3 2 20 1 1 1 1 4 1 1 1 5 3 2 1 6 1 4'
a.=3
Do i=1 By 1 while al<>''
Parse Var al a.i al
End
pix=calc(194)
Do e=1 To length(pi)
If substr(pix,e,1)<>substr(pi,e,1) Then Leave
End
Numeric Digits 50
Say 'PI ='||(pi+0)||'...'
Say 'PIX='||(pix+0)||'...'
Say (e-1) 'correct digits'
Exit
 
Get_Coeffs: procedure Expose a b a.
Parse Arg n
a=a.n
b=1
Return
 
Calc: procedure Expose a b a.
Parse Arg n
Temp=0
do ni = n to 1 by -1
Call Get_Coeffs ni
Temp = B/(A + Temp)
end
call Get_Coeffs 0
return (A + Temp)

[edit] 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')
Output:
$ ruby cfrac.rb                                                              
1.41421356237309504880168872420969807856967187537695
2.71828182845904523536028747135266249775724709369996
3.1415926536

[edit] Scala

Works with: Scala version 2.9.1

Note that Scala-BigDecimal provides a precision of 34 digits. Therefore we take a limitation of 32 digits to avoiding rounding problems.

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()
}
}
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:

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()
}
}
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

[edit] 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.

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
Output:
1.4142135623730951
2.7182818284590455
3.1415926373965735

[edit] XPL0

The number of iterations (N) needed to get the 13 digits of accuracy was determined by experiment.

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);
]
Output:
16
1.414213562372820
1.414213562373100
13
2.718281828459380
2.718281828459050
10000
3.141592653589540
3.141592653589790

[edit] zkl

fcn cf(fa,fb,a0){fcn(fa,fb,a0,n){
a0 + [n..1,-1].reduce('wrap(p,n){fb(n)/(fa(n)+p)},0.0) }
.fp(fa,fb,a0) }

cf creates a function that calculates the continued fraction from the bottom up. The new function takes a single parameter, n, which is used to calculate the nth term.

sqrt2:=cf((2.0).noop,(1.0).noop,1.0);
sqrt2(200) : "%.20e".fmt(_).println();
nap:=cf((0.0).create,fcn(n){(n==1) and 1.0 or (n-1).toFloat()},2.0);
println(nap(15) - (1.0).e);
pi:=cf((6.0).noop,fcn(n){n=2*n-1; (n*n).toFloat()},3.0);
println(pi(1000) - (1.0).pi);

(1.0).create(n) --> n, (1.0).noop(n) --> 1.0

Output:
1.41421356237309514547e+00
1.33227e-15
-2.49251e-10
Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox