[分享] 計算 e 到小數點下一億位

[Schottky (Schottky)]

[9/17 12:15pm 更新: 增加位數範圍檢查, 修正 terms() 問題。]

這是為了完成一個小時候的夢想... 計算 e 到小數點下一億位...
                   ∞
所使用的級數是 e = Σ(1/k!) = 1/0! + 1/1! + 1/2! + 1/3! + 1/4! + .....
                  k=0
這個級數看似簡單, 其實很優秀, 計算不難, 且收斂速度非常快。

程式碼: http://codepad.org/qtZaidZq

這次一樣用二分法, 把 N 項的級數切成一半, 左右分別算完再合併結果,
要注意的是需要計算的總項數和分數通分加法數量和迴圈法其實是一樣多的,
那為什麼遞迴 Divide & Conquer 會比迴圈法快這麼多呢?
是因為上次我打馬虎眼被 Feis 抓包的大數運算複雜度並不是 O(1) ...
二分法切割到最後都是小小的運算, 兩項兩項小數字合併, 速度就很快,
迴圈法跑沒多久, 加總的數字開始變得很大, 就舉步維艱了...

不要被 Divide & Conquer 這個名詞騙了, 重點在 Conquer 之後的 merge,
只要能找到有效率的 merge 公式就成功了... 幸運地, 在這邊是有的,
把計算結果用分數表示 (分子&分母分別用整數表示) 再通分即可合併。

千萬不要全部轉成浮點數計算, 每一次運算都用小數點以下一億位的超高精度浮點
運算會非常慢... 整數的話只有最後的幾次超級合併才有這麼大的數量級...
而且化成分數最後只進行一次超級除法, 可以把浮點誤差降到最低。

一開始需要估計計算到第幾項才能達到我們要的精度。
因為不知道怎麼找出 Stirling's approximation 的反函數, 我以二分逼近法求解。

在 64-bit Linux, 2.4GHz CPU 上面大約需要兩分半鐘, 所以我寫了進度條,
其中有一分鐘是將最後 120MB 左右的計算結果轉換成十進位字串, 以及寫到輸出檔。
如果一開始就使用十進位(十億進位)來儲存和運算, 最後的轉換過程就可以省下來。

我在 Linux 和 MinGW 上都 compile 過, 32-bit 環境也可以跑, 沒問題的。
輸入位數最大可以到 2^64-1, 問題是 RAM 不夠, 一億(10^8)位就要 2GB RAM 了,
估計要計算十億位需要 20GB RAM, 各位如果手邊有暴力機台可以測試看看 :)

// URL: http://codepad.org/qtZaidZq
// -------- >8 -------- >8 -------- >8 -------- >8 -------- >8 --------
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <gmp.h>
#include <math.h>

#define FILENAME    "e.txt"

#define PI  3.141592653589793238462643L
#define E   2.718281828459045235360287L

unsigned long long int count=0, total=0, old_p=0;

long double logfactorial(unsigned long long int n)
{
    return (logl(2.0*PI)/2.0+logl(n)/2.0+n*logl(n/E))/logl(10.0);
}

unsigned long long int terms(unsigned long long int digits)
{
    long double f;
    unsigned long long int  n, upper, lower;

    for (n=1;;n<<=1) {
        f = logfactorial(n);
        if (f>digits) break;
    }
    upper = n;
    lower = 0;
    while (upper-lower > 1) {
        n = (upper+lower)/2;
        f = logfactorial(n);
        if (f>digits) {
            upper = n;
        } else {
            lower = n;
        }
    }
    n = upper;

    return n;
}

void write_e(mpf_t e, unsigned long long int digits)
{
    char *str;
    long int    exp;
    int c;
    FILE *fp;

    fp = fopen(FILENAME, "w+");
    if (fp == NULL) { return; }

    str = malloc(digits+2);

    mpf_get_str(str, &exp, 10, digits+1, e);

    fprintf(fp, "  e = ");
    for (c=0; c<digits; c++) {
        if (c!=1) {
            if (c%50 == 1) {
                fprintf(fp, "\t");
            }
        }
        fprintf(fp, "%c", str[c]);
        if (c==0) {
            fprintf(fp, ".");
        } else if (c%1000 == 0) {
            fprintf(fp, " << %d\n", c);
        } else if (c%500 == 0) {
            fprintf(fp, " <\n");
        } else if (c%50 == 0) {
            fprintf(fp, "\n");
        } else if (c%5 == 0) {
            fprintf(fp, " ");
        }
    }

    if (c%50 != 1) {
        fprintf(fp, "\n");
    }

    free(str);

    fclose(fp);
}

void progress(void)
{
    int p, g, c;

    p = (int)floorf(1000.0f*count/total+0.5f);
    if (p > old_p) {
        old_p = p;
        g = (int)floorf(72.5f*count/total+0.5f);
        for (c=0; c<72; c++) {
            if (c<g) printf("H");
            else printf("-");
        }

        printf(" %5.1f%%\r", p/10.0);
        fflush(stdout);
    }

    if (count == total) printf("\n");
}

void calc(mpz_t a, mpz_t b, mpz_t p, mpz_t q)
{
    mpz_t   m, p2, q1;

    mpz_inits(m, p2, q1, NULL);

    mpz_sub(m, b, a);
    if (mpz_cmp_ui(m, 0) <= 0) {
        if (mpz_cmp_ui(a, 0) == 0) {
            mpz_set_ui(p, 1);
        } else {
            mpz_set_ui(p, 0);
        }
        mpz_set_ui(q, 1);
    } else if (mpz_cmp_ui(m, 1) == 0) {
        if (mpz_cmp_ui(a, 0) == 0) {
            mpz_set_ui(p, 2);
            mpz_set_ui(q, 1);
        } else {
            mpz_set_ui(p, 1);
            mpz_set(q, b);
        }
    } else {
        mpz_add(m, a, b);
        mpz_fdiv_q_ui(m, m, 2);
        calc(a, m, p, q1);
        calc(m, b, p2, q);

        // Merge
        mpz_mul(p, p, q);
        mpz_add(p, p, p2);

        mpz_mul(q, q, q1);
    }

    //gmp_printf("P(%d,%d)=%Zd, Q(%d,%d)=%Zd\n", a, b, p, a, b, q);
    count++;
    progress();

    mpz_clears(m, p2, q1, NULL);
}

void calc_e(unsigned long long int digits)
{
    mpz_t   z, p, q, n;
    mpf_t   pf, qf;
    unsigned long long int d, n_terms, precision;

    if (digits < ULLONG_MAX) {
        d = digits + 1;
    } else {
        d = ULLONG_MAX;
    }
    n_terms = terms(d);
    total = n_terms*2-1;
    precision = ceill(d*(logl(10.0L)/logl(2.0L)))+1;

    printf("d = %llu, n = %llu, precision = %llu\n", d, n_terms, precision);
    mpf_set_default_prec(precision);
    precision = mpf_get_default_prec();
    printf("Real precision = %llu\n", precision);

    mpz_inits(z, p, q, n, NULL);
    mpf_inits(pf, qf, NULL);

    mpz_set_ui(z, 0);
    mpz_set_ui(n, n_terms);
    calc(z, n, p, q);

    mpf_set_z(pf, p);
    mpf_set_z(qf, q);
    mpf_div(pf, pf, qf);

    mpz_clears(z, p, q, n, NULL);

    printf("Calculation done. Writing result to file...\n");
    write_e(pf, d);

    mpf_clears(pf, qf, NULL);
}

int main(int argc, char *argv[])
{
    unsigned long long int digits = 100000000ULL;
    unsigned long long int prec_max;

    // I don't know how GNU MP sets it's maximum precision, it's black magic.
    if (ULONG_MAX == 4294967295UL) {
        prec_max = 1292913975UL;    // 32-bit
    } else {
        prec_max = 5553023288523357111ULL;  // 64-bit
    }

    if (argc == 2) {
        digits = strtoull(argv[1], NULL, 10);
        // 1292913985 = MAX_PRECISION*log10(2)-1 , MAX_PRECISION = 2^32-1
        if (digits == 0 || digits > prec_max) {
            fprintf(stderr, "Invalid parameter, "
                "digits must range from 1 to %llu\n", prec_max);
            return -1;
        }
    } else if (argc > 2) {
        fprintf(stderr, "Usage: %s #digits\n", argv[0]);
        return -1;
    }

    calc_e(digits);

    return 0;
}

--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 220.137.9.131
推 LPH66:先推再說
推 damody:我以為是硬幹 原來是mpfr。
→ damody:跑 terms 函數 數字就爆了= =
感謝提醒, terms() 的確是疏忽了。

這裡用的是 GNU MP 不是 MPFR, 原因是 mpf_ 系列運算稍微快一點點,
而且都已經用了 gmp, 我不想再多寫一份 mpfr 安裝說明...
事實上 gmp/mpfr 兩個版本我都有寫, 所以差不多可以另外寫一篇評比了 =.="

「硬幹」的意思是? 如果是指 implement 大數運算, 那整個程式碼會膨脹很多,
但想要突破 RAM 的限制把 HD 拿來當計算紙, gmp 是辦不到的, 只能自己寫。
目前世界紀錄 10^13 位數就是這樣做的, HD 當 RAM 用, RAM 當 cache 用。

推 f814030:太酷了！
推 lc85301:好猛喔LOL
※ 編輯: Schottky        來自: 220.137.7.154        (09/17 12:33)
推 MOONRAKER:Great
推 damody:了解 詳細解說~