ガウス＝ルジャンドル法に基づく円周率の計算プログラム(C言語版)

Gauss–Legendre.c

#define N 100

#include<stdio.h>
#include<math.h>

int main()
{
  int i=0;
  double a[N+1];
  double b[N+1];
  double t[N+1];
  double p[N+1];
  double pi = 0.00 ;
  a[0] = 1.0;
  b[0] = 1.0/sqrt(2);
  t[0] = 1.0/4;
  p[0] = 1;
  for(i=1;i<N+1;i++)
  {
      a[i] = (a[i-1]+b[i-1])/2;
      b[i] = sqrt(a[i-1]*b[i-1]);
      t[i] = t[i-1] - (p[i-1]*(a[i-1]-a[i])*(a[i-1]-a[i]));
      p[i] = 2 * p[i-1];
  }
        pi = ((a[N]+b[N])*(a[N]+b[N]))/(4*t[N]);
          
          printf("%30f\n",pi);
          return 0;          
        }