% Archimedes1.m
% Compute pi by finding the perimeters of inscribing 
% regular polygons of sides 2^(n+1), n = 1, 2, ...
% iteratively for large n.
% This version uses a recursion relation for tn.
% Algoritm involves one lesss square root and fewer
% number of multiplications and divisions.
% The savings in computational time is also measured.
% However It still suffers from bad cancellations.
% K. Ming Leung, 01/28/04

clear; format long; numberLoops = 10000;
% version 0
tic;
for loop = 1:numberLoops
twoN1=4; 
n1 = 1;
P(n1)= sqrt(8);
for n=1:30
    n1 = n+1;
    P(n1) = twoN1*sqrt(2*(1-sqrt(1-(P(n)/twoN1)^2)));
    %disp([n1 P(n1)]);
    twoN1=twoN1*2;    
end;
end;
toc;

tic;
for loop = 1:numberLoops
half = 0.5; 
n1 = 1;
%twoN1 = 4;
T(n1)= half;
%disp([n1 P(n1)]);
for n=1:30
    n1 = n+1;
    T(n1) = half*(1-sqrt(1-T(n)));
    %twoN1=twoN1*2;
    %disp([n1 twoN1*sqrt(T(n1))]);
end;
end;
toc
% Iterating tn instead of pn is about 40% faster!
