/* 円周率計算プログラム */
/* マチンの公式         */
/* pi = 16 * arctan( 1/5 ) - 4 * arctan( 1/239 ) */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define MAX_T 150000L
#define CARRY 10000000L
/* １ロングワード＝１０進７桁格納  */
/* 最大１０５００００桁（１００万桁）まで計算可能（でもやると死ぬ） */

long StartTime, EndTime;
long w1[MAX_T], w2[MAX_T], p[MAX_T], s[MAX_T];

void devine( ary, dat, stp, len )
  long      *ary, dat, stp, len;
{
register long i, r;

r = 0;
for( i = stp; i < len; i++ )
  {
  register long w;
  register double m;
  
  m = CARRY * (double)r + ary[i];
  w = ary[i] = m / dat;
  r = m - (double)dat * w;
  }
}

void main( argc, argv )
  int argc;
  char *argv[];
{
register long w1p, w2p, n, l;
register short f;

if( argc == 1 )
  {
  printf("PAI [length]\n");
  exit(2);
  }

l = atol( argv[1] ) / 7L + 3L;
if( (l > MAX_T) || (l < 4L) )
  {
  printf("PAI parameter error\n");
  exit(1);
  }

time( &StartTime );

for( n = 0L; n < l; n++ )
  {
  w1[n] = w2[n] = p[n] = 0L;
  }

f = 0;
w1[1] = 80L;  /* 16 * 5   */
w2[1] = 956L; /*  4 * 239 */
w1p = 1L;
w2p = 1L;
for( n = 1L; ; n += 2L )
  {
  register long i, r;
  
  devine( &w1[0], 25L, w1p, l );
  if( f == 0 )
    {
    devine( &w2[0], 57121L, w2p, l );
    if( w2[w2p] == 0L )
      w2p++;
    if( w2p >= l )
      f = 1;
    }
  r = 0;
  for( i = l - 1L; i >= w1p; i-- )
    {
    s[i] = w1[i] - w2[i] - r;
    if( s[i] < 0L )
      {
      s[i] += CARRY;
      r = 1L;
      }
    else
      r = 0L;
    }
  devine( &s[0], n, w1p, l );
  if( (n & 2L) == 2L )
    {
    register long j;
    
    for( j = l - 1L; j >= 1L; j-- )
      {
      p[j] -= s[j];
      if( p[j] < 0L )
        {
        p[j] += CARRY;
        --p[j-1];
        }
      }
    }else{
    register long j;

    for( j = l - 1L; j >= 1L; j-- )
      {
      p[j] += s[j];
      if( p[j] >= CARRY )
        {
        p[j] -= CARRY;
        ++p[j-1];
        }
      }
    }
  if( w1[w1p] == 0L )
    w1p++;
  if( w1p >= l )
    break;
  }

time( &EndTime );

printf("PAI= %ld.\n", p[1] );
f = 100;
for( n = 2L; n < l; n++ )
  {
  switch( f )
    {
    case 7:
      printf("%07lu\n", p[n] );
      break;
    case 6:
      printf("%06lu\n%01lu", p[n]/10L, p[n] % 10L );
      break;
    case 5:
      printf("%05lu\n%02lu", p[n]/100L, p[n] % 100L );
      break;
    case 4:
      printf("%04lu\n%03lu", p[n]/1000L, p[n] % 1000L );
      break;
    case 3:
      printf("%03lu\n%04lu", p[n]/10000L, p[n] % 10000L );
      break;
    case 2:
      printf("%02lu\n%05lu", p[n]/100000L, p[n] % 100000L );
      break;
    case 1:
      printf("%01lu\n%06lu", p[n]/1000000L, p[n] % 1000000L );
      break;
    default:
      printf("%07lu", p[n] );
      break;
    }
  f -= 7;
  if( f <= 0 )
    f += 100;
  }
printf("\n計算時間=%lu sec\n", EndTime-StartTime );
}
/* end of file */