quick ugly port of bellard.org/pi/pi1.c

["Ruud H.G. van Tol" <rvtol@isolution.nl>]

Just for fun, I started a port of <https://bellard.org/pi/pi1.c>.
To my surprise it appeared to work at the first try.
Plenty of bug fixes and optimizations still to do!


DIVN( t, a, v, vinc, kq, kqinc )= {
  kq+= kqinc;
  if(kq >= a
  , while( kq >= a, kq-= a);
    if( kq == 0
    , while( !(t % a)
      , t/= a;
        v+= vinc;
      )
    )
  );
  [t, v, kq];
}

main(n)= {
  my( av, a, vmax, N, num, den, k, kq1, kq2, kq3, kq4, t, v, s, i, sum=0);

  N= floor((n + 20) * log(10) / log(13.5));

  a= 1;
  while( a <= 3*N
  , a= nextprime(a+1);
    vmax= floor(log(3*N) / log(a));
    if( a == 2
    , vmax+= (N - n);
      vmax>0 || next;
    );

    av= a^vmax;
    s= 0;
    den= 1;
    kq1= 0;
    kq2= -1;
    kq3= -3;
    kq4= -2;

    if (a == 2
    , num= 1;
      v= -n
    , num= 2^n % av;
      v= 0
    );

    for( k=1, N
    , t= 2 * k;
      [t, v, kq1]= DIVN(t, a, v, -1, kq1, 2);
      num= num * t % av;

      t= 2 * k - 1;
      [t, v, kq2]= DIVN(t, a, v, -1, kq2, 2);
      num= num * t % av;

      t= 3 * (3 * k - 1);
      [t, v, kq3]= DIVN(t, a, v, 1, kq3, 9);
      den= den * t % av;

      t= (3 * k - 2);
      [t, v, kq4]= DIVN(t, a, v, 1, kq4, 3);
      if(a != 2, t*= 2, v++);
      den= den * t % av;

      if( v > 0
      , if( a != 2
        , if( av%2, t= 1/den % av)
        , t= 1/den % av
        );
        t= t * num % av;
        for( i=v, vmax-1, t= t * a % av);
        t= t * (25 * k - 3) % av;
        s+= t;
        if( s >= av, s-= av );
      )
    );

    t= 5^(n - 1) % av;
    s= s * t % av;
    sum= (sum + s / av) % 1.0;
  );

  printf("Decimal digits of pi at position %d: %09d\n"
  , n
  , floor(sum * 1e9)
  );
  return(0);
}



Example runs:

? main(1)
Decimal digits of pi at position 1: 141592653
cpu time = 7 ms, real time = 7 ms.

? main(761)
Decimal digits of pi at position 761: 499999983
cpu time = 1,487 ms, real time = 1,498 ms.

-- Ruud