function value = agm_pi ( d )

%*****************************************************************************80
%
%% agm_pi(): Arithmetic-geometric mean for pi.
%
%  Discussion:
%
%    The Brent-Salamin algorithm is used to compute PI to D decimal digits.
%
%    a0 = 1
%    b0 = 1 / sqrt ( 2 )
%    t0 = 1 / 4
%    p0 = 1
%
%    repeat:
%      an+1 = ( an + bn ) / 2
%      bn+1 = sqrt ( an bn )
%      tn+1 = tn - pn * ( an - an+1 )^2
%      pn+1 = 2 *pn
%
%    pi is approximated by ( an+1 + bn+1 )^2 / 4 / tn+1
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    27 February 2012
%
%  Author:
%
%    Cleve Moler
%
%  Reference:
%
%    Cleve Moler,
%    Cleve's Corner, "Computing Pi",
%    https://www.mathworks.com/company/newsletters/news_notes/2011/ 
%
%  Input:
%
%    integer D, the number of decimal digits desired.
%
%  Output:
%
%    symbolic VALUE, the value of pi to D digits.
% 

%
%  VPA items will use D digits.
%
  digits ( d );
%
%  Initial conditions:
%    a = 1
%    b = 1 / sqrt(2)
%
  a = vpa ( 1, d );
  b = 1 / sqrt ( vpa ( 2, d ) );
  s = 1 / vpa ( 4, d );
  p = 1;
%
%  Accuracy increases quadratically.
%  Estimate number of iterations needed.
%
  n = ceil ( log2 ( d ) );

  for k = 1 : n
    c = ( a + b ) / 2;
    b = sqrt ( a * b );
    s = s - p * ( c - a )^2;
    p = 2 * p;
    a = c;
  end

  value = a^2 / s;

  return
end
 
