C++ - ネイピア数（自然対数の底）e 計算！

以前、コンピュータで大きな桁数を計算する概念・アルゴリズムを紹介しました。

• C++ - 多桁計算！

今回は、ネイピア数（自然対数の底） \(e\) を多桁計算するアルゴリズムについてです。

当然、プログラミン言語そのものが保有している関数は使用しません。級数展開して計算します。

0. 前提条件

• Linux Mint 14 Nadia (64bit) での作業を想定。
• g++ (Ubuntu/Linaro 4.7.2-2ubuntu1) 4.7.2

また、当方の環境で扱える int 型の範囲は -2,147,483,648 〜 2,147,483,647 であることから、

• 多桁計算で使用する1つの配列のサイズは8桁としている。

1. 計算概要

ネイピア数 \(e\) をマクローリン展開（テイラー展開において \(x=1\) としたもの）すると、 \(\displaystyle e = 1 + \frac{1}{1!} + \frac{1}{2!} + \frac{1}{3!} + \cdots + \frac{1}{n!} + \cdots\) となる。また、小数点以下 1,000 桁まで計算する場合、 \(449! < 10^{1000} < 450!\) であるので、余裕を見て、分母に階乗のある項を第 451 項まで計算することにする。

また、ネイピア数は何種類か定義方法がありますが、ここでは説明しません。

2. C++ ソース作成

例として、以下のようにソースを作成した。
計算する桁数は、小数点以下 1,000 桁としている。
ソース内の L の値を変更し、対応する計算項数 N も変更すれば、任意の桁数を計算可能。

File: calc_napier.cpp

/*********************************************
 * ネイピア数 e (自然対数の底)計算
 *********************************************/
#include <iostream>  // for cout
#include <stdio.h>   // for printf()

#define L  1000           // 算出桁数
#define L1 ((L / 8) + 1)  // 配列サイズ
#define L2 (L1 + 1)       // 配列サイズ + 1
#define N  451            // 計算項数

using namespace std;

/*
 * 計算クラス
 */
class Calc
{
    // 各種変数
    int k, i;           // LOOP インデックス
    int carry, borrow;  // 繰り上がり、借り
    long w;             // 被乗数、被除数ワーク
    long r;             // 剰余ワーク

    public:
        // 計算
        void calc();
        // ロング + ロング
        void ladd(int *, int *, int *);
        // ロング / ショート
        void ldiv(int *, int, int *);
        // 結果出力
        void display(int *);
};

/*
 * 計算
 */
void Calc::calc()
{
    // 配列宣言
    int s[L2 + 2], a[L2 + 2];

    // 配列初期化
    s[0] = a[0] = 1;
    for (k = 1; k <= L2; k++)
        s[k] = a[k] = 0;

    // 計算
    for (k = 1; k <= N; k++) {
        ldiv(a, k, a);
        ladd(s, a, s);
    }

    // 結果出力
    display(s);
}

/*
 * ロング + ロング
 */
void Calc::ladd(int a[], int b[], int c[])
{
    carry = 0;
    for (i = L2; i >=0; i--) {
        c[i] = a[i] + b[i] + carry;
        if (c[i] < 100000000) {
            carry = 0;
        } else {
            c[i] -= 100000000;
            carry = 1;
        }
    }
}

/*
 * ロング / ショート
 */
void Calc::ldiv(int d[], int e, int f[])
{
    r = 0;
    for (i = 0; i < L2; i++) {
        w = d[i];
        f[i] = (w + r) / e;
        r = ((w + r) % e) * 100000000;
    }
}

/*
 * 結果出力
 */
void Calc::display(int s[])
{
    printf("%7d. ", s[0]);
    for (i = 1; i < L1; i++)
        printf("%08d ", s[i]);
    printf("\n");
}

/*
 * メイン処理
 */
int main()
{
    try
    {
        // 計算クラスインスタンス化
        Calc objCalc;

        // ネイピア数計算
        objCalc.calc();
    }
    catch (...) {
        cout << "例外発生！" << endl;
        return -1;
    }

    // 正常終了
    return 0;
}

• Gist - C++ source code to compute Napier’s constant.

3. C++ ソースコンパイル

$ g++ calc_napier.cpp -o calc_napier

何も出力されなければ成功です。

4. 実行

$ ./calc_napier
      2. 71828182 84590452 35360287 47135266 24977572 47093699 95957496 
69676277 24076630 35354759 45713821 78525166 42742746 63919320 03059921 
81741359 66290435 72900334 29526059 56307381 32328627 94349076 32338298 
80753195 25101901 15738341 87930702 15408914 99348841 67509244 76146066 
80822648 00168477 41185374 23454424 37107539 07774499 20695517 02761838 
60626133 13845830 00752044 93382656 02976067 37113200 70932870 91274437 
47047230 69697720 93101416 92836819 02551510 86574637 72111252 38978442 
50569536 96770785 44996996 79468644 54905987 93163688 92300987 93127736 
17821542 49992295 76351482 20826989 51936680 33182528 86939849 64651058 
20939239 82948879 33203625 09443117 30123819 70684161 40397019 83767932 
06832823 76464804 29531180 23287825 09819455 81530175 67173613 32069811 
25099618 18815930 41690351 59888851 93458072 73866738 58942287 92284998 
92086805 82574927 96104841 98444363 46324496 84875602 33624827 04197862 
32090021 60990235 30436994 18491463 14093431 73814364 05462531 52096183 
69088870 70167683 96424378 14059271 45635490 61303107 20851038 37505101 
15747704 17189861 06873969 65521267 15468895 70350354 

当方の非力なマシンでも１秒かかりませんでした。
自分の環境と相談して、計算する桁数を大きくしてみるのもよいでしょう。

但し、結果が何万桁以上になる場合は、多桁（多倍長）乗算に上記の方法ではなく「Karatsuba法」や「Toom-Cook法」や「FFT（高速フーリエ変換）」を使うのが一般的なようだ。

---

やはり、「多桁（多倍長）計算は計算機に限る」というを実感できます。

※ちなみに最近の当方の C++ アルゴリズムについての記事は、古い C によるアルゴリズムに関する書物を参考に C++ に移植した形態となっています。

以上。