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

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

On 2024-12-14 15:17, Ruud H.G. van Tol wrote:
> Just for fun, I started a port of <https://bellard.org/pi/pi1.c>.
> [...]
> Plenty of bug fixes and optimizations still to do!

Slightly optimized it:

pi_dec(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;
  until( 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;
      if( (kq1+=2) >= a
      , until( kq1 < a, kq1-= a );
        kq1 || until( t%a, t\= a; v--)
      );
      num= num * t % av;

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

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

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

      if( v > 0
      , if( a == 2
        , t= 1/den % av
        , if( av%2, 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);
}


Run examples:

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


? pi_dec(761)
Decimal digits of pi at position 761: 499999983
cpu time = 793 ms, real time = 802 ms.


? pi_dec(2000)
Decimal digits of pi at position 2000: 994657640
cpu time = 4,539 ms, real time = 4,578 ms.


To check, see https://oeis.org/A000796/b000796.txt
which has 20k entries, but its index is +1.

-- Ruud