function [a,b]=agm(a0,b0,N)

% AGM
% a(n+1) = (a(n)+b(n))/2;
% b(n+1) = sqrt(a(n)*b(n));
% To connect to K(sin alpha), let a0 = 1, b0 = cos alpha
% Should give pi/[2 M(1,cos alpha)]
% Also, AGM(a,b) = pi/4*(a+b)/ellipticK( (a-b)/(a+b) )

% Examples:
% agm(1,2,5);
% agm(1,sqrt(2),4);  % Gauss' constant = 1/agm(1,sqrt(2))
% agm(24,6,5);
% agm(1,1/sqrt(2),4);  % Use for pi approximation

a(1) = a0;
b(1) = b0;
for n=1:N-1
    a(n+1) = (a(n)+b(n))/2;
    b(n+1) = sqrt(a(n)*b(n));
end

xpy=a(N)+b(N);
xmy=a(N)-b(N);
M=pi/4*xpy/ellipticK(xmy/xpy); % relate AGM to 

disp(['a = ', num2str(a(N)), ' b = ', num2str(b(N)), ' M = ', num2str(M)])
disp(['a-b =', num2str(a(N)-b(N)), '  a-M =', num2str(a(N)-M)])


% Gauss-Legendre Algorithm for pi
% Typing
% agm(1,1/sqrt(2),5);
% gives the results:
%    a = 0.84721 b = 0.84721 M = 0.84721
%    a-b =1.1102e-16  a-M =0
%    pi = 3.1416    error = -8.8818e-16
% To displsay more digits, enter the following:
% format long

t(1) = 1/4;
p(1) = 1;
for n=1:N-1
    t(n+1) = t(n)-p(n)*(a(n)-a(n+1))^2;
    p(n+1) = 2*p(n);
end
P = xpy^2/4/t(N);
% disp(['pi = ', num2str(P), '    error = ', num2str(pi-P)])


