% Archimedes0.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 the perimeters.
% K. Ming Leung, 01/28/04

clear; format long;

twoN1=4; 
n1 = 1;
P(n1)= sqrt(8);
disp([n1 P(n1)]);
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

% After nZero = log2(3*sqrt(4/eps)) iterations, P should go
% to zero because of rounding and cancellation.
nZero = log2(3*sqrt(4/eps)),

