// Mat01169, script para Metodo das Medias para Pi 
function d=digse(ra,r)
 if (r==ra) | (abs(r)<%eps)
    d=16.0;   // maximo que pode ser obtido em IEEE 754 //
   // com precisao dupla (por enquanto...) //
 else
    d=max(0,-log10( 2*abs( (ra-r)/r ))); // digse //
 end
endfunction; // funcao que calcula os digses //

// Calculo de Pi pelo Metodo das Medias (Arquimedes c. 200 AC) //
nmax=20;
ai=3*sqrt(3); bi=3*sqrt(3)/2;  // inicializacao de a1 e b1 //
printf("%i %20.16f %20.16f %10.4f \n",1,ai,bi,0); 
for n=2:nmax
  aa=ai; ba=bi;
  ai = 2/(1/aa + 1/ba);  // media harmonica  //
  bi = sqrt(ai*ba);      // media geometrica //
  printf("%i %20.16f %20.16f %10.4f \n",n,ai,bi,digse(aa,ai)); 
end
// criado Ago 2004, modificado Ago 2011 //



