fork download
copy 
  1. /*
  2.  * Computation of the n'th decimal digit of \pi with very little memory.
  3.  * Written by Fabrice Bellard on January 8, 1997.
  4.  * 
  5.  * We use a slightly modified version of the method described by Simon
  6.  * Plouffe in "On the Computation of the n'th decimal digit of various
  7.  * transcendental numbers" (November 1996). We have modified the algorithm
  8.  * to get a running time of O(n^2) instead of O(n^3log(n)^3).
  9.  * 
  10.  * This program uses mostly integer arithmetic. It may be slow on some
  11.  * hardwares where integer multiplications and divisons must be done
  12.  * by software. We have supposed that 'int' has a size of 32 bits. If
  13.  * your compiler supports 'long long' integers of 64 bits, you may use
  14.  * the integer version of 'mul_mod' (see HAS_LONG_LONG).  
  15.  */
  16.  
  17. #include <stdlib.h>
  18. #include <stdio.h>
  19. #include <math.h>
  20.  
  21. /* uncomment the following line to use 'long long' integers */
  22. /* #define HAS_LONG_LONG */
  23.  
  24. #ifdef HAS_LONG_LONG
  25. #define mul_mod(a,b,m) (( (long long) (a) * (long long) (b) ) % (m))
  26. #else
  27. #define mul_mod(a,b,m) fmod( (double) a * (double) b, m)
  28. #endif
  29.  
  30. /* return the inverse of x mod y */
  31. int inv_mod(int x, int y)
  32. {
  33.     int q, u, v, a, c, t;
  34.  
  35.     u = x;
  36.     v = y;
  37.     c = 1;
  38.     a = 0;
  39.     do {
  40.     q = v / u;
  41.  
  42.     t = c;
  43.     c = a - q * c;
  44.     a = t;
  45.  
  46.     t = u;
  47.     u = v - q * u;
  48.     v = t;
  49.     } while (u != 0);
  50.     a = a % y;
  51.     if (a < 0)
  52.     a = y + a;
  53.     return a;
  54. }
  55.  
  56. /* return (a^b) mod m */
  57. int pow_mod(int a, int b, int m)
  58. {
  59.     int r, aa;
  60.  
  61.     r = 1;
  62.     aa = a;
  63.     while (1) {
  64.     if (b & 1)
  65.         r = mul_mod(r, aa, m);
  66.     b = b >> 1;
  67.     if (b == 0)
  68.         break;
  69.     aa = mul_mod(aa, aa, m);
  70.     }
  71.     return r;
  72. }
  73.  
  74. /* return true if n is prime */
  75. int is_prime(int n)
  76. {
  77.     int r, i;
  78.     if ((n % 2) == 0)
  79.     return 0;
  80.  
  81.     r = (int) (sqrt(n));
  82.     for (i = 3; i <= r; i += 2)
  83.     if ((n % i) == 0)
  84.         return 0;
  85.     return 1;
  86. }
  87.  
  88. /* return the prime number immediatly after n */
  89. int next_prime(int n)
  90. {
  91.     do {
  92.     n++;
  93.     } while (!is_prime(n));
  94.     return n;
  95. }
  96.  
  97. int main(int argc, char *argv[])
  98. {
  99.     int av, a, vmax, N, n, num, den, k, kq, kq2, t, v, s, i;
  100.     double sum;
  101.  
  102.     if (argc < 2 || (n = atoi(argv[1])) <= 0) {
  103.     printf("This program computes the n'th decimal digit of \\pi\n"
  104.            "usage: pi n , where n is the digit you want\n");
  105.     exit(1);
  106.     }
  107.  
  108.     N = (int) ((n + 20) * log(10) / log(2));
  109.  
  110.     sum = 0;
  111.  
  112.     for (a = 3; a <= (2 * N); a = next_prime(a)) {
  113.  
  114.     vmax = (int) (log(2 * N) / log(a));
  115.     av = 1;
  116.     for (i = 0; i < vmax; i++)
  117.         av = av * a;
  118.  
  119.     s = 0;
  120.     num = 1;
  121.     den = 1;
  122.     v = 0;
  123.     kq = 1;
  124.     kq2 = 1;
  125.  
  126.     for (k = 1; k <= N; k++) {
  127.  
  128.         t = k;
  129.         if (kq >= a) {
  130.         do {
  131.             t = t / a;
  132.             v--;
  133.         } while ((t % a) == 0);
  134.         kq = 0;
  135.         }
  136.         kq++;
  137.         num = mul_mod(num, t, av);
  138.  
  139.         t = (2 * k - 1);
  140.         if (kq2 >= a) {
  141.         if (kq2 == a) {
  142.             do {
  143.             t = t / a;
  144.             v++;
  145.             } while ((t % a) == 0);
  146.         }
  147.         kq2 -= a;
  148.         }
  149.         den = mul_mod(den, t, av);
  150.         kq2 += 2;
  151.  
  152.         if (v > 0) {
  153.         t = inv_mod(den, av);
  154.         t = mul_mod(t, num, av);
  155.         t = mul_mod(t, k, av);
  156.         for (i = v; i < vmax; i++)
  157.             t = mul_mod(t, a, av);
  158.         s += t;
  159.         if (s >= av)
  160.             s -= av;
  161.         }
  162.  
  163.     }
  164.  
  165.     t = pow_mod(10, n - 1, av);
  166.     s = mul_mod(s, t, av);
  167.     sum = fmod(sum + (double) s / (double) av, 1.0);
  168.     }
  169.     printf("Decimal digits of pi at position %d: %09d\n", n,
  170.        (int) (sum * 1e9));
  171.     return 0;
  172. }
Runtime error #stdin #stdout 0s 3456KB
comments (?)
stdin
copy 
pi 5 5 5 5 5
stdout
copy 
This program computes the n'th decimal digit of \pi
usage: pi n , where n is the digit you want
https://ideone.com/J4Nqs1
language:
C++ (gcc 8.3)
created:
10 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!