C 言語サンプルプログラム (円周率 10000 桁)

-- 目次 --

  1. 前書き
  2. arctan(1/m)
  3. 円周率 1000 桁
  4. 円周率 10000 桁
  5. πの古典的公式

1. 前書き

古代から円周率は色々な方法で計算され、 現在では特別な方法を使うと、 円周率を 10 億桁以上も計算することができます。 但し、これはスーパコンピュータでないと無理な方法です。 ここでは、パソコンで円周率を計算するときによく使用されるマチンの公式から 単純に円周率を求めるプログラムを考えます。

なお長倍精度演算に高速フーリエ変換と呼ばれる手法を取り入れると、 非常に速く計算することができますが、 円周率 1 万桁程度であれば (大したことがない数値ですから)、 特殊な方法を取り入れなくても瞬く間にパソコンで計算することができます。 何も知らない人であれば多分息を飲むことになると思います。

最初にマチンの公式を思い出しておこう。

arctan1
m
=
1
m
-  1 
3 m3
+  1 
5 m5
-  1 
7 m7
+ . . .   (m 自然数)
 π =
4 ( 4 arctan 1
5
- arctan  1 
239
)

2. arctan(1/m)

arctan(1/m) を求めることができるプログラムを作り、 バグを取るために arctan(1/5) の計算をして、Windows の電卓で計算した arctan(1/5) の値と比較することにします。 そのため最初は桁数をあまり大きくしません。

atan.c

#include <stdio.h>

#define M 8     // 4*(M-1)=小数点以下の桁数 (10 進表示)
#define N 7     // 4*(N-1)=小数点以下の桁数 (16 進表示)
                // M-1≒1.204(N-1)        (M, N の関係)

void Add(unsigned a[], unsigned b[]);
void Sub(unsigned a[], unsigned b[]);
void Div(unsigned a[], unsigned d);
void Mul(unsigned a[], unsigned d);
void Dup(unsigned a[], unsigned b[]);
void ATAN(unsigned x[], int m);
int if_zero(unsigned a[]);
void Display(unsigned a[]);
void to_decimal(unsigned a[], unsigned w[]);



void main(){
    unsigned x[N];

    ATAN(x, 5);
    Display(x);
}
// arctan(1/m) を x に格納

void ATAN(unsigned x[], int m){
    int i, e;
    unsigned m2;
    unsigned y[N], z[N];
    for(i = 0; i < N; i++){     // 配列 x,y の初期化
        x[i] = 0;
        y[i] = 0;
    }
    y[0] = 1;

    Div(y, m);
    Dup(y, x);
    m2 = m*m;
    i = 1;
    do {
        Div(y, m2);
        Dup(y, z);
        Div(z, 2*i+1);
        if( (i & 1) == 0 ){     // i が偶数のとき
            Add(x, z);
        } else {                // i が奇数のとき
            Sub(x, z);
        }
        i++;
        e = if_zero(z);         // z が 0 でないと -1
    } while (e == -1);          // z が 0 でない限り継続
    
}

//足算,  a + b を a に代入
void Add(unsigned a[], unsigned b[]){
    int i,j;
    unsigned x;
    for(i = 0; i < N; i++){
        x = a[i] + b[i];
        if(x <= 0xffff){            // 桁上がりしないとき
            a[i] = x;
        } else {                    // 桁上がりするとき
            a[i] = x & 0xffff;      // x の下位 2 バイト
            j = i-1;                // 以下桁上がりの操作
            while(a[j] == 0xffff){  // 1 加えると更に桁が上がるとき
                a[j] = 0;
                j--;
            }
            a[j]++;
        }
    }
}

//引き算,  a - b を a  に代入 (但し a>b)
void Sub(unsigned a[], unsigned b[]){
    int i,j;
    for(i = 0; i < N; i++){
        if (a[i] >= b[i]){                  // そのまま引けるとき
            a[i] = a[i] - b[i];
        } else {                            // 上の桁から借りるとき
            a[i] = 0x10000 + a[i] - b[i];   // 1 を上から借りる
            j = i-1;
            while(a[j] == 0){               // 上の桁が 0 であれば
                a[j] = 0xffff;              // 更に上から借りる
                j--;
            }
            a[j]--;
        }
    }

}

//割り算,  a/d を a に代入 (d <= 0xffff, 16 進 4桁)
// 上の桁から順番に d で割り算。
void Div(unsigned a[], unsigned d){
    int i;
    unsigned x, q, res;
    res = 0;
    for(i = 0; i < N; i++){
        res = res << 16;        // 上の桁の余りを 16 進で4 桁ずらして
        x = a[i] + res;         // 加える
        q = x/d;                // d でわる
        a[i] = q;               // 商
        res = x - q*d;          // 余りは下の桁の割り算に加える
    }
}

//かけ算, a*d を a に代入 (d <= 0xffff, 16 進 4 桁)
// 下の桁から順番に d をかける。
void Mul(unsigned a[], unsigned d){
    int i; 
    unsigned x, q;
    q = 0;
    for(i = N-1; i >= 0; i--){
        x = a[i]*d + q;         // 各桁の d 倍に桁上がりの部分を加える
        a[i] = x & 0xffff;      // x の下位 2 バイト (16 進 4 桁)
        q = x >> 16;            // x の上位 2 バイト (16 進 4 桁)
    }                           // 桁上がりの部分
}

//a を b にコピー
void Dup(unsigned a[], unsigned b[]){
    int i;
    for(i = 0; i < N; i++){
        b[i] = a[i];
    }
}

// a = 0 のとき 0 を返し、そうでないとき -1 を返す関数
int if_zero(unsigned a[]){
    int i, j;
    j = 0;
    for(i = 0; i < N; i++){
        if (a[i] != 0) {            // 一箇所でも 0 でないとき
            j = -1;                 // j の値を -1 にする
            break;
        }
    }
    return j;
}

//  10 進数に変換して表示
void Display(unsigned a[]){
    int i,j;
    unsigned w[M];
    to_decimal(a, w);               // 16 進数 a を 10 進数 w に変換
    printf("%4.1u.",w[0]);          // 小数点以上の表示
    for(i = 1; i < M; i++){
        printf("%4.4u ", w[i]);     // 4 桁ごと表示
        j = i%12;                   // j は i を 12 で割った余り
        if(j==0) printf("\n     "); // 12 個表示したとき改行
    }
    printf("\n");                   // 一番最後に改行
}

//  10 進数への変換
void to_decimal(unsigned a[], unsigned w[]){
    int i;
    unsigned b[N];
    for(i = 0; i < N; i++){
        b[i] = a[i];
    }
    w[0] = b[0];
    b[0] = 0;
    for(i = 1; i < M; i++){
        Mul(b, 10000);              // b を 10000 倍して
        w[i] = b[0];                // 小数点以上を w[i] に代入
        b[0] = 0;                   // 小数点以上を 0 にする
    }
}

プログラムの簡単な説明

大まかな説明をすると、大きさ N の配列は 16 進小数をあらわし、 大きさ M の配列は 10 進小数をあらわす。各関数の機能は

ATAN(x, m) actan(1/m) を配列 x に格納する関数
Add(a, b) 配列 a (16 進小数とみなす) と配列 b の和を配列 a に格納する関数
Sub(a, b) 配列 a から配列 b を引き、配列 a に格納する関数。配列 a が 配列 b より大きいと仮定している。
Div(a, b) 配列 a を 16 進数 d (0≦d<216) で割って、商を配列 a に格納する関数。
Mul(a, d) 配列 a に 16 進数 d (0≦d<216) をかけて、積を配列 a に格納する関数。
Dup(a, b) 配列 a を配列 b にコピーする関数
if_zero(a) 配列 a が 0 であれば 0 を返し、そうでなければ -1 を返す関数。
Display(a) 16 進数 a を 10 進数に変換して表示する。この関数から以下の関数が呼び出される。
to_decimal(a, w) 16 進小数 a を 10 進小数に変換して w に格納する関数。

ATAN(x,m) の関数内部では x,y,z は次のように変化する

y : 
1
m
,  1 
m3
,  1 
m5
, . . .
z : 
1
m
,  1 
3m3
,  1 
5m5
, . . .
z : 
1
m
, 1
m
- 1 
3m3
, 1
m
- 1 
3m3
+  1 
5m5
, . . .

16 進数の表示

プログラム中で 16 進数を使う場合には先頭に 0x (ゼロ、エックス) をつける。例えば


0xffff
 

は 164 - 1 (= 216 - 1) のことです。

ビットごとの演算子

ビットとは 2 進数の一桁のことで、ビットごとの演算子には例えば次のようなものがある。

演算子意味
&AND
|OR
<<左シフト
>>右シフト

シフト命令は説明ずみですから、AND(&) と OR(|) の説明をします。例えば 4 桁の 2 進数は 次のようにして 4 個の要素 {a1, a2, a3, a4} からなる集合 X の部分集合と見なせます。


(0001)2 <--> {a1},   (0010)2 <--> {a2},   (1010)2 <--> {a4, a2}
 

このとき X の部分集合としての共通部分 (AND) と和集合 (OR) に対応するのが ビットごとの AND と OR です。具体的には

(0011)2 & (0001)2 --> (0001)2
(0011)2 | (0001)2 --> (0011)2
(0111)2 & (0101)2 --> (0101)2
(0111)2 & (0101)2 --> (0111)2

あるいは 1 桁だけに注目すれば

1 & 1 --> 1    1 | 1 --> 1
1 & 0 --> 0    1 | 0 --> 1
0 & 0 --> 0    0 | 0 --> 0

従って、次のような使い方ができます (ビット演算は符号なし整数の場合のみに使用を 限定する方がよい)。

  1. a, b が符号なし整数のとき

    
    b = a & 1;
    

    とすれば, a が奇数のときに b が 1 になり、a が偶数のときに b は 0 になる。

  2. a, b が符号なし整数のとき

    
    b = a & 0xffff;
    

    とすれば a を 2 進表示したときの下 16 桁 (下 2 バイト) が b に代入されます。

論理演算は次のことを実現します。わざと右側にビット演算を対応させます。

&& 真 --> 真   1 & 1 --> 1
&& 偽 --> 偽   1 & 0 --> 0
&& 偽 --> 偽   0 & 0 --> 0
|| 真 --> 真   1 | 1 --> 1
|| 偽 --> 真   1 | 0 --> 1
|| 偽 --> 偽   0 | 0 --> 0

みてわかるように、真は値 1, 偽は値 0 です。実は関係演算子は、 真のときは 1 を返し、偽のときは 0 を返しています。 これは次のようなプログラムで検証することができます。


#include <stdio.h>

void main()
{
     int a,b,c;

    printf("整数を入力してください :");
    scanf("%d", &a);

    printf("整数を入力してください :");
    scanf("%d", &b);

    c = (a==b);
    printf("%d\n", c);
}

同じ整数を入力したときに 1 が表示され、 異なる整数を入力したときに 0 が表示される。

関係演算子は真の場合には値 1 を返すが、0 以外の値は C 言語では実はすべて真として定義されている。 これは C 言語の標準仕様 (ANSI C) である。 (こうなっている理由は機械語 (アセンブラー) では 0 に等しい場合とそうでない場合にわけて 制御することが多いため、C 言語でもそのように定義する方が都合がよいのであろう。)

プログラムのチェック

Windows NT/95/98 には電卓 (アクセサリーの中) がついているから


ATAN(x, 5);

の 5 の部分をいろいろ取り替えて、電卓の計算結果と一致することを確認する。 Windows の電卓に関して注意する点は

下の絵は Windows 95 の電卓で、クリックする場所を示している。 他の Windows のバージョンでは若干位置が違っている。

3. 円周率 1000 桁

arctan(1/m) のプログラムを使って、円周率を 1000 桁求めることにしましょう。 誤差の評価をすれば一番下の数桁が違っているだけであることが分かりますから、 誤差のことは気にしないことにします。 ともかくも円周率を求めるために atan.c のメイン関数を次のように変更します。


void main(){
    unsigned x[N], y[N];

    ATAN(x, 5);
    Mul(x, 4);

    ATAN(y, 239);

    Sub(x, y);
    Mul(x, 4);

    Display(x);
}

更に 1000 桁求めるために


M = 251 (= 1000/4 + 1)
 

とおき、16 進小数10 進小数と 16 進小数の関係 から


N≒ 250÷1.204 + 1≒208.6
 

となるから N = 209 とおこう。従って


#define M 251   // 4*(M-1)=小数点以下の桁数 (10 進)
#define N 209   // 4*(N-1)=小数点以下の桁数 (16 進)
                // M-1≒1.204 (N-1)

としよう。以上によって円周率 1000 桁が求まる。 出力が以下の結果とほぼ同じことになることを確認しよう

円周率 1000 桁 (岩波数学辞典から)

   3.1415 9265 3589 7932 3846 2643 3832 7950 2884 1971 6939 9375 
     1058 2097 4944 5923 0781 6406 2862 0899 8628 0348 2534 2117 
     0679 8214 8086 5132 8230 6647 0938 4460 9550 5822 3172 5359 
     4081 2848 1117 4502 8410 2701 9385 2110 5559 6446 2294 8954 
     9303 8196 4428 8109 7566 5933 4461 2847 5648 2337 8678 3165 
     2712 0190 9145 6485 6692 3460 3486 1045 4326 6482 1339 3607 
     2602 4914 1273 7245 8700 6606 3155 8817 4881 5209 2096 2829 
     2540 9171 5364 3678 9259 0360 0113 3053 0548 8204 6652 1384 
     1469 5194 1511 6094 3305 7270 3657 5959 1953 0921 8611 7381 
     9326 1179 3105 1185 4807 4462 3799 6274 9567 3518 8575 2724 
     8912 2793 8183 0119 4912 9833 6733 6244 0656 6430 8602 1394 
     9463 9522 4737 1907 0217 9860 9437 0277 0539 2171 7629 3176 
     7523 8467 4818 4676 6940 5132 0005 6812 7145 2635 6082 7785 
     7713 4275 7789 6091 7363 7178 7214 6844 0901 2249 5343 0146 
     5495 8537 1050 7922 7968 9258 9235 4201 9956 1121 2902 1960 
     8640 3441 8159 8136 2977 4771 3099 6051 8707 2113 4999 9998 
     3729 7804 9951 0597 3173 2816 0963 1859 5024 4594 5534 6908 
     3026 4252 2308 2533 4468 5035 2619 3118 8171 0100 0313 7838 
     7528 8658 7533 2083 8142 0617 1776 6914 7303 5982 5349 0428 
     7554 6873 1159 5628 6388 2353 7875 9375 1957 7818 5778 0532 
     1712 2680 6613 0019 2787 6611 1959 0921 6420 1989 

4. 円周率 10000 桁

以上のプログラムで 10 進数で小数点以下何桁まで arctan(1/m) を計算できるかを考えてみよう。

今 10 進数で小数点以下 L 桁まで求めるとすると arctan(1/m) の和で

 1 
2i+1
×  1  
m2i+1
< 1 
10L
(*)

となる最初の i までの和を加えなければならない。ところがプログラム中の Div(a, b) は長倍精度の a を 2 バイトの b で割ることしかできません。つまり (*) を満たす最小の i に対して


2i+1 < 216 (**)

となる必要がある。(*) より


log10(2i+1) + (2i+1) log10 m > L
 

(**) より


log10(2i+1) < log10 216 = log10 65536 < 5
 

となるから、L はほぼ


L < 216 log10 m
 

の範囲である。m = 5 とすれば


L < 216 log10 5 = 45807.6982
 

従って、以上のプログラムでは arctan(1/5) は 45000 桁程度しか計算できないことになる。 当然この範囲であれば


#define M 251     // 4*(M-1)=小数点以下の桁数 (10 進)
#define N 209     // 4*(N-1)=小数点以下の桁数 (16 進)
                  // M-1≒1.204(N-1)

の M, N の値を変えてもプログラムは正常に動作するはずである。

印刷のことなど、リダイレクト

10000 桁も出力すれば、画面が流れて、読み取れなくなる。 円周率のプログラムが pi.c ですでにビルドしてあれば、実行ファイル pi.exe ができていて、通常はDOS 窓から


pi

と打ち込めば、pi.exe が実行されるが、もしも


pi >pi.txt

と打ち込むと、出力結果がすべてファイル pi.txt に出力される。 このファイルはテキストエディターでもワープロでも読み出すことができ、 印刷することが可能である。

printf標準出力 に出力するためのコマンドで、 標準出力は通常はディスプレーであるが ">" によって標準出力を 変更することが可能である。これを出力のリダイレクト (redirect = 向きを変えること) という。入力に関しても同様である。

計算チェック

円周率を 10000 桁出力しても、結果が正しいかどうか分からない。 このような場合には別の方法で円周率を求めて、その比較をすることから計算チェックを すればよい。 円周率の公式にはマチンの公式以外にも色々な公式がある。そこで上で与えたプログラムの メイン関数を少々手直しをして、別の公式から円周率を計算して比較すれば 計算チェックが可能となる。

5. π の古典的公式

π の古典的公式には次のようなものがある。

π=
16 arctan1
5
- 4 arctan  1 
239
(Machin)
 =
24 arctan 1
8
+ 8 arctan  1 
57
+ 4 arctan  1 
239
(Stoemer)
 =
48 arctan  1 
18
+ 32 arctan  1 
57
- 20 arctan  1 
239
(Gauss)
 =
32 arctan  1 
10
- 4 arctan  1 
239
- 16 arctan  1 
515
(Klingenstiema)

これまでの方法では割り算に関しては多倍長の小数を 4 桁以内の 16 進数 (= 2 バイト 整数) で割る方法しか与えていないため、求められる桁数には制限がある。 この点から最も効率のよいのが Gauss の公式で、このとき


L < 216 log10 18 = 82265. . . .
 

であるから 8 万桁程度は求められることになる。