円周率とネイピア数

たまには"技研"らしいエントリでも書いてみましょっか。
以前フェルマー＝ワイルズの定理について触れましたが（フェルマー＝ワイルズの定理 - しまうま技研）、証明に重要な役割を果たした谷山＝志村＝ヴェイユ予想（現在では谷山＝志村定理）について、ウィキペディアでは以下のように記述されています。

谷山・志村定理（または谷山・志村予想）の意味は、「楕円曲線論」と「保型形式論」という異なる二つの分野で用いられる特殊な概念が同種のものである、ということだ。これは、予想が提出された当時では、とても衝撃的なことだった。なぜなら、この二つの分野のそれぞれの概念は、まったく別の概念だと思われていたからだ。

二つの分野の別の概念が同種のものだとすれば、そこには何らかの深遠な真実がひそんでいることになる。それゆえ、この予想は、通常の定理のように一つの分野だけの問題ではなくて、数学における広範な真実を告げる重大な問題だと理解された。

谷山–志村予想 - Wikipedia

フェルマー予想は分かりやすいけれど、数学的には大した意味はないと考えられていました。
実際、天才数学者ガウスは「フェルマーの定理のような、証明も反証もできないような式はいくらでも書き出せる」なんていっていますしね。（cf. フェルマーの最終定理）
一方の谷山＝志村予想は（提起直後はともかく）数学上の重大な問題として捉えられていて、それがフェルマー予想を解く重要な鍵になったのは何となく皮肉めいたものを感じます。

さて、谷山＝志村定理の解説にある「この二つの分野のそれぞれの概念は、まったく別の概念だと思われていた」という一文を読んで、私は真っ先にオイラーの等式を連想しました。

• オイラーの公式 - Wikipedia

それぞれ全く異なる定義から導き出される定数が、虚数単位を介して結び付いている不思議な式。　以下のように書かれることもあります。

こう書くと、加法の単位元 0 と乗法の単位元 1 が現れて、より一層ミステリアスさが引き立つとか。

さてこの円周率ですが、自分で計算したことがある人はどのくらいいますかね？
今なら Superπ などをダウンロードしてきてチャッチャと計算しちゃうのでしょうが、私のようなオッサンが"マイコン"に触れ始めた頃は、みんな自分で工夫してプログラミングしてたんですよね。
んで、Superπのソースを探してみたんだけど、ちょっと調べてみた限りではバイナリしか見当たりませんでした。
なので、昔カシオのFX-602Pで組んだプログラムをベースに、Cで書き直して見ましたよ。

/*************************************/
/*                                   */
/*  **  PI calculation program  **   */
/*                                   */
/* Copyright (C)2007 Zebra Research  */
/*   http://d.hatena.ne.jp/k87p561/  */
/*                                   */
/*************************************/

#include <stdio.h>
#include <malloc.h>
#include <memory.h>

#define LIMIT 1000  /* 計算桁数の上限（デフォルト値） */

/***************/
/* 多桁加算処理 */
/***************/
int add( signed char *dst, signed char *src, int figure )
{
    int carry = 0;
    int result = 0;
    int i;

    for ( i = figure - 1; i >= 0; i-- ) {
        dst[ i ] += ( src[ i ] + carry );
        if ( dst[ i ] > 9 ) {
            dst[ i ] -= 10;
            carry = 1;
        }
        else {
            carry = 0;
        }
        result += src[ i ];
    }
    return result;
}

/***************/
/* 多桁減算処理 */
/***************/
int sub( signed char *dst, signed char *src, int figure )
{
    int borrow = 0;
    int result = 0;
    int i;

    for ( i = figure - 1; i >= 0; i-- ) {
        dst[ i ] -= ( src[ i ] + borrow );
        if ( dst[ i ] < 0 ) {
            dst[ i ] += 10;
            borrow = 1;
        }
        else {
            borrow = 0;
        }
        result += src[ i ];
    }
    return result;
}

/***************/
/* 多桁除算処理 */
/***************/
void div( signed char *dst, int n, int figure )
{
    int mod = 0;
    int tmp;
    int i;

    for ( i = 0; i < figure; i++ ) {
        mod *= 10;
        tmp = mod + dst[ i ];
        if ( tmp >= n ) {
            dst[ i ] = tmp / n;
            mod = tmp % n;
        }
        else {
            mod += dst[ i ];
            dst[ i ] = 0;
        }
    }
}

/*********************/
/* 多桁逆正接計算処理 */
/*********************/
void Arctan( signed char *dst, int k, int x, int figure )
{
    signed char *x2;
    signed char *x3;
    int n = 1;
    int sign = 1;
    int result;

    /* 一時計算領域の確保 */
    x2 = ( signed char * )malloc( figure );
    x3 = ( signed char * )malloc( figure );

    memset( dst,  0, figure );
    memset( x2, 0, figure );
    memset( x3, 0, figure );

    x2[ 0 ] = k;
    div( x2, x, figure );

    for ( ; ; ) {
        memcpy( x3, x2, figure ); /* backup */

        div( x2, n, figure );

        if ( sign > 0 ) {
            result = add( dst, x2, figure );
        }
        else {
            result = sub( dst, x2, figure );
        }

        /*******************************/
        /* ループ終了条件判定           */
        /*                             */
        /* 加算・減算する値が 0 なら終了 */
        /*******************************/
        if ( result == 0)
            break;

        n += 2;
        sign *= ( -1 );

        memcpy( x2, x3, figure ); /* restore */

        /* x^2 での除算                            */
        /*   2回に分けて除算することで精度を確保する */
        div( x2, x, figure );
        div( x2, x, figure );
    }

    /* 一時計算領域の解放 */
    free( x2 );
    free( x3 );
}

/*********************/
/* calculation of PI */
/*********************/
void main( int argc, char *argv[] )
{
    signed char *pi;
    signed char *x1;
    int figure;
    int i;

    /* 引数で指定された桁数をチェック＆セット */
    if ( argc == 2 ) {
        if ( sscanf( argv[ 1 ], "%d", &figure ) != 1 ) {
            figure = LIMIT;
        }
        else if ( figure <= 0 ) {
            figure = LIMIT;
        }
    }
    else {
        figure = LIMIT;
    }

    /* 計算領域の確保 */
    pi = ( signed char * )malloc( figure );
    x1 = ( signed char * )malloc( figure );

    /* 計算処理 */
    Arctan( pi, 16, 5, figure );  /* pi = 16 * arctan( 1/5 ) */

    Arctan( x1, 4, 239, figure ); /* x1 = 4 * arctan( 1/239 ) */
    sub( pi, x1, figure );        /* pi = 16 * arctan( 1/5 ) - 4 * arctan( 1/239 ) */

    /* 結果の表示 */
    printf( "PI=%1d.\n\n", pi[ 0 ] ); /* 1桁目の表示 */
    for ( i = 1; i < figure; i++ ) {
        printf( "%1d", pi[ i ] );
        if ( i % 10 == 0 ) {
            if ( i % 50 == 0 ) {
                if ( i % 1000 == 0 ) { /* 1000桁毎に空行挿入（ブロック化） */
                    printf( "\n\n" );
                }
                else {
                    printf( "\n" ); /* 50桁毎に改行 */
                }
            }
            else {
                printf( " " ); /* 10桁毎の区切り */
            }
        }
    }
    printf( "\n" );

    /* 計算領域の解放 */
    free( pi );
    free( x1 );
}

これはマーチンの公式（マチンの公式とも）を使って円周率を求めるプログラム。　一応、Microsoft C/C++ で動作確認しました。　第1引数で指定した桁数分、円周率を計算してくれます。（引数を省略したり、無効な値が指定された時は、デフォルトで1000桁が適用されます）
ちなみに、マーチンの公式とは、以下のようなもの。

プログラムでは式を下記のように変形して使用しています。

収束の早いマーチンの公式を用いることや多桁計算のアイディアは、はるか昔、工学社のI/Oに掲載された記事から拝借しました。
またこのプログラムを元にして、割と簡単にネイピア数を求めるプログラムに書き直せます。
Cコンパイラをお持ちの方は、暇つぶしにどうぞ。
なお、int型は 32bit であることを期待した作りになっていますので、16bit の処理系をお使いの方は一部 long にした方が良いかも知れません。

最後にお約束の使用許諾条件を列挙しておきます。

• このプログラムの転載・改変はご自由にどうぞ。　いちいち私（作者）に断りを入れる必要もありません。
• このプログラムを使用したことによる損害等については、作者は一切関知致しません。
• このプログラムの処理やメモリ使用効率はあまり良いものではありません。　Superπなどと比べて貰っても困ります。

あ、そうそう、計算結果の検算なんかは、こちらのサイトをご覧になって下さいな。

• πの部屋！