/* 円周率計算プログラム */
/* 高野の公式         */
/* pi = 48*arctan(1/49)+128*arctan(1/57)-20*arctan(1/239)+48*arctan(1/110443) */

#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], w3[MAX_T], w4[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, w3p, w4p, n, l;
register short f;

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

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

time( &StartTime );

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

w1[1] = 2352L;    /*  48 * 49     */
w2[1] = 7296L;    /* 128 * 57     */
w3[1] = 4780L;    /*  20 * 239    */
w4[1] = 5301264L; /*  48 * 110443 */
w1p = 0L;
w2p = 1L;
w3p = 1L;
w4p = 1L;
for( n = 1L; ; n += 2L )
  {
  register long i, r;
  
  devine( &w1[0], 2401L, w1p, l );
  if( w2p < l )
    {
    devine( &w2[0], 3249L, w2p, l );
    if( w2[w2p] == 0L )
      w2p++;
    }
  if( w3p < l )
    {
    devine( &w3[0], 57121L, w3p, l );
    if( w3[w3p] == 0L )
      w3p++;
    }
  if( w4p < l )
    {
    devine( &w4[0], 110443L, w4p, l );
    devine( &w4[0], 110443L, w4p, l );
    if( w4[w4p] == 0L )
      w4p++;
    }
  if( w4p < l )
    {
    r = 0L;
    for( i = l - 1L; i >= w1p; i-- )
      {
      register long m;

      m = w1[i] + w2[i] - w3[i] + w4[i] + r;
      if( m < 0L )
        {
        r = -1L;
        s[i] = m + CARRY;
        }
      else
        {
        for( r = 0L; m >= CARRY; r++ )
          m -= CARRY;
        s[i] = m;
        }
      }
    }
  else if ( w3p < l )
    {
    r = 0L;
    for( i = l - 1L; i >= w1p; i-- )
      {
      register long m;

      m = w1[i] + w2[i] - w3[i] + r;
      if( m < 0L )
        {
        r = -1L;
        s[i] = m + CARRY;
        }
      else
        {
        for( r = 0L; m >= CARRY; r++ )
          m -= CARRY;
        s[i] = m;
        }
      }
    }
  else if ( w2p < l )
    {
    r = 0L;
    for( i = l - 1L; i >= w1p; i-- )
      {
      register long m;

      m = w1[i] + w2[i] + r;
      if( m >= CARRY )
        {
        r = 1L;
        s[i] = m - CARRY;
        }
      else
        {
        r = 0L;
        s[i] = m;
        }
      }
    }
  else
    {
    for( i = l - 1L; i >= w1p; i-- )
      s[i] = w1[i];
    }
  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 */