ninedigitsofpiat.cs

中文名：Windows Forms 程序设计 英文名：Windows Forms Programming in c# 作者： Chris Sells 翻译： 荣耀 蒋贤哲 出版社：人民

// NineDigitsOfPiAt.cs

/*
 * Computation of the n'th decimal digit of pi with very little memory.
 * Written by Fabrice Bellard on January 8, 1997.
 * Ported to C# by Chris Sells on May 5, 2002.
 * 
 * We use a slightly modified version of the method described by Simon
 * Plouffe in "On the Computation of the n'th decimal digit of various
 * transcendental numbers" (November 1996). We have modified the algorithm
 * to get a running time of O(n^2) instead of O(n^3log(n)^3).
 * 
 * This program uses mostly integer arithmetic. It may be slow on some
 * hardwares where integer multiplications and divisons must be done
 * by software.
 */

using System;

namespace MySecondApp {

  public class NineDigitsOfPi {

    public static int StartingAt(int n) {
      int av = 0;
      int vmax = 0;
      int N = (int)((n + 20) * Math.Log(10) / Math.Log(2));
      int num = 0;
      int den = 0;
      int kq = 0;
      int kq2 = 0;
      int t = 0;
      int v = 0;
      int s = 0;
      double sum = 0.0;

      for( int a = 3; a <= (2 * N); a = next_prime(a) ) {
        vmax = (int)(Math.Log(2 * N) / Math.Log(a));
        av = 1;

        for( int i = 0; i < vmax; ++i ) av = av * a;

        s = 0;
        num = 1;
        den = 1;
        v = 0;
        kq = 1;
        kq2 = 1;

        for( int k = 1; k <= N; ++k ) {
          t = k;
          if( kq >= a ) {
            do {
              t = t / a;
              --v;
            }
            while( (t % a) == 0 );

            kq = 0;
          }

          ++kq;
          num = mul_mod(num, t, av);

          t = (2 * k - 1);
          if( kq2 >= a ) {
            if( kq2 == a ) {
              do {
                t = t / a;
                ++v;
              }
              while( (t % a) == 0 );
            }

            kq2 -= a;
          }

          den = mul_mod(den, t, av);
          kq2 += 2;

          if( v > 0 ) {
            t = inv_mod(den, av);
            t = mul_mod(t, num, av);
            t = mul_mod(t, k, av);
            for( int i = v; i < vmax; ++i ) t = mul_mod(t, a, av);
            s += t;
            if( s >= av ) s -= av;
          }
        }

        t = pow_mod(10, n - 1, av);
        s = mul_mod(s, t, av);
        sum = (sum + (double)s / (double)av) % 1.0;
      }

      return (int)(sum * 1e9);
    }

    private static int mul_mod(long a, long b, int m) {
      return (int)((a * b) % m);
    }

    // return the inverse of x mod y
    private static int inv_mod(int x, int y) {
      int q = 0;
      int u = x;
      int v = y;
      int a = 0;
      int c = 1;
      int t = 0;

      do {
        q = v / u;

        t = c;
        c = a - q * c;
        a = t;

        t = u;
        u = v - q * u;
        v = t;
      }
      while( u != 0 );

      a = a % y;
      if( a < 0 ) a = y + a;

      return a;
    }

    // return (a^b) mod m
    private static int pow_mod(int a, int b, int m) {
      int r = 1;
      int aa = a;

      while( true ) {
        if( (b & 0x01) != 0 ) r = mul_mod(r, aa, m);
        b = b >> 1;
        if( b == 0 ) break;
        aa = mul_mod(aa, aa, m);
      }

      return r;
    }

    // return true if n is prime
    private static bool is_prime(int n) {
      if( (n % 2) == 0 ) return false;

      int r = (int)(Math.Sqrt(n));
      for( int i = 3; i <= r; i += 2 ) {
        if( (n % i) == 0 ) return false;
      }

      return true;
    }

    // return the prime number immediately after n
    private static int next_prime(int n) {
      do {
        n++;
      }
      while( !is_prime(n) );

      return n;
    }
  }
}