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

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

Ported it a bit deeper, used Mods:

pi_dec(n)= {
  my( a=1, sum=0, N );

  n>0||n=1;
  N= floor((n + 20) * log(10) / log(13.5));

  until( a > 3*N
  , my( kq1=0, kq2=-1, kq3=-3, kq4=-2, av, den, num, s, v, vmax );

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

    av= a^vmax;

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

    for( k=1, N
    , my(t);

      t= 2 * k;
      if( (kq1+=2) >= a
      , (kq1 % a) || until( t%a, t\= a; v-- )
      );
      num*= t;

      t= 2 * k - 1;
      if( (kq2+=2) >= a
      , (kq2 % a) || until( t%a, t\= a; v-- )
      );
      num*= t;

      t= 9 * k - 3;
      if( (kq3+=9) >= a
      , (kq3 % a) || until( t%a, t\= a; v++ )
      );
      den*= t;

      t= 3 * k - 2;
      if( (kq4+=3) >= a
      , (kq4 % a) || until( t%a, t\= a; v++ )
      );
      if( a == 2, v++, t*= 2);
      den*= t;

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

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

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

  return(0);
}

-- Ruud