#Header Begins--Do not remove##################################################
#
#	Maple code to compute pi a multitude of different ways.
#
#	(c) Eugene Fiume, University of Toronto, 1995 (modified 1996).
#
#	You are permitted to use this code for noncommercial, educational
#	purposes only.  This header must not be removed from this file.
#
#Header Ends--Do not remove####################################################

#
#	Leibnitz' rule to get for pi/4.
#
L := proc (n)
	local i;
	sum( (-1)^i / (2*i+1), i=0..n );
end:

#
#	Sequence to compute Pi by Leibnitz' rule.
#
P1 := proc (n)
	local i;
	4*sum( (-1)^i / (2*i+1), i=0..n );
end:
#
#	Sequence to compute Pi by Euler's rule.
#
P2 := proc (n)
	local i;
	sqrt(6*sum( 1 / i^2, i=1..n ));
end:


#
#	Compute pi by perimeter of regular 2^n-gon inscribing a circle of
#	unit diameter.  Thus P(2) is 2*sqrt(2).  (Archimedes' rule.)
#
P3 := proc (n)
	option remember;
	if   (n<2) then 0
	elif (n=2) then 2*sqrt(2)
	else 2^(n-1)*sqrt(2*(1-sqrt(1-(P3(n-1)/2^(n-1))^2)));
	fi;
end:


#
#	Compute pi by perimeter of regular n-gon inscribing a circle of
#	unit diameter.  This is the nonrecursive form that uses Pi to
#	compute it!
#
P4 := proc (n)
	local i;
	n * sqrt( 1/2*(1-cos(2*Pi/n)) );
end:


#
#	The recursive form of P4 to avoid use of Pi.  See the notes.
#
P4r  := proc (n)
	local i;
	2^n*sqrt(1/2*(1-C(n-1)));
end:

C    := proc (n)
	option remember;
	if   (n<2) then 0
	elif (n=2) then sqrt(2)/2
	else sqrt( 1/2*(C(n-1)+1) );
	fi;
end:


#
#	This evaluates any of the above schemes to a floating point
#	representation.  It is useful to observe how sensitive the
#	above schemes are to floating point errors.
#
Pe := proc (P,n)
	evalf(P(n));
end:


#
#	Compute a sequence of values of one of the sequences above suitable
#	for plotting or computing errors.
#
compP := proc(P,n)
	local i;

	[seq( [i, P(i)], i=4..n)];
end:


#
#	Plot a sequence of approximations to pi using one of the above.
#
plotP := proc(P,n)
	plot(compP(P,n));
end:


#
#	The routines below compute various forms of error between the
#	sequence approximation and Pi.
#
errP := proc(P,n)
	local v,p;
	p   := compP(P,n);
	v   := Pi;
	plot([ [i, v-p[i][2]]$i=1..nops(p)]);
end:

#
#	Absolute error for approximation of procedure P to value v.
#
errPabs := proc(P,v,n)
	local p;
	p   := compP(P,n);
	plot([ [i, abs(v-p[i][2])]$i=1..nops(p)]);
end:

#
#	Relative error for approximation of procedure P to value v.
#
errPrel := proc(P,v,n)
	local p;
	p   := compP(P,n);
	plot([ [i, abs(v-p[i][2])/v]$i=1..nops(p)]);
end:

#
#	Do composite absolute error for some of the above routines.
#	(This is mainly just an example of how to superimpose several
#	plots into one display.)  You can also use the "display"
#	procedure to do this.
#
errPabsAll := proc(N)
	local p1,p2,p3,p4,p5, e1,e2,e3,e4,e5, n;
	option remember;
	p1  := compP(P4r,N);
	p4  := compP(P4,N);
	e1  := [ [i, abs(Pi-p1[i][2])]$i=1..nops(p1)];
	e4  := [ [i, abs(Pi-p4[i][2])]$i=1..nops(p4)];

	plot({e1,e4},n=0..nops(p1));
end:

#
#	Sample executions:
#
#	P1(10);    # (or P2, or P3, or ...)
#
#	Pe(P1,100);
#
#	Compare the following:
#
#	evalf(P4r(10));
#	P3(P4r,10);
#
#	Why are these results different?
#
#	plotP(P4r,10);
#	errPrel(P4r,15);   # explain the results!!