計算誤差とはそういうもの

> options(digits=20)
> (x <- sum(1/factorial(0:18)))
[1] 2.718281828459045
> x == exp(1)
[1] TRUE
> (y <- 1+sum(1/factorial(1:18)))
[1] 2.718281828459046
> y == exp(1)
[1] FALSE

ちなみに，計算精度の範囲内では，

> (x <- sum(1/factorial(0:17)))
[1] 2.718281828459045
> x == exp(1)
[1] TRUE

で，十分

ちなみに，奥村先生のアルゴリズム入門にある C プログラムは長い。
数学の教科書に書かれていることを確認してみようかなあと思う人も，C で確かめようとすると大変だけど，R だと超簡単。

/***********************************************************
    e.c -- 自然対数の底
***********************************************************/

long double ee(void)
{
    int n;
    long double e, a, prev;

    e = 0;  a = 1;  n = 1;
    do {
        prev = e;  e += a;  a /= n;  n++;
    } while (e != prev);
    return e;
}

#include <stdio.h>
#include <stdlib.h>
#include <float.h>

int main()
{
    printf("e = %.*Lg¥n", LDBL_DIG, ee());
    return EXIT_SUCCESS;
}