为什么我的C程序用于计算质量差的欧拉常数？

#include <stdio.h>

long factorial(int n) {
    long result = 1;
    for (int i = 1; i <= n; ++i)
    result *= i;
    return result;
}

int main ()
{
    double n=0;
    int i;
    for (i=0; i<=32; i++) {
        n=(1.0/factorial(i))+n;
    }
    printf("%.32f\n", n);
}

#1 楼

如果要实现这样的事情，则应首先了解如何完成这些事情。我希望这听起来不太苛刻。为了更好地解释，这是我的代码变体，与C使用expl(1)计算的结果进行了比较：

#include <stdio.h>
#include <math.h>

int main ()
{
    long double n = 0, f = 1;
    int i;
    for (i = 28; i >= 1; i--) {
        f *= i;  // f = 28*27*...*i = 28! / (i-1)!
        n += f;  // n = 28 + 28*27 + ... + 28! / (i-1)!
    }  // n = 28! * (1/0! + 1/1! + ... + 1/28!), f = 28!
    n /= f;
    printf("%.64llf\n", n);
    printf("%.64llf\n", expl(1));
    printf("%llg\n", n - expl(1));
    printf("%d\n", n == expl(1));
}

输出： >

2.7182818284590452354281681079939403389289509505033493041992187500
2.7182818284590452354281681079939403389289509505033493041992187500
0
1

您的代码有两个重要更改：


此代码不计算1，1 * 2，1 * 2 * 3,。 ..是O（n ^ 2），但是一次计算1 * 2 * 3 * ...（这是O（n））。

它从较小的数字开始。让我们暂时假设您的阶乘是正确的（见下文）。当您计算

1/1 + 1/2 + 1/6 + ... + 1/20！

并尝试将其添加到1/21！时，您正在添加

1/21！ = 1/51090942171709440000 = 2E-20，

到2.something，对结果没有影响（双精度保留约16个有效数字）。这种影响称为下溢。

但是，您是否从这些数字开始，即，如果您计算出1/32！+1/31！+ ...都会产生一定的影响。


注意32！需要118位，即15个字节。您的long类型占用的空间不大（标准要求至少32位；容纳的位数不可能超过64）。如果将printf("%ld\n", factorial(i));添加到您的主for -loop中，则会看到阶乘：

1
2
6
24
120
720
5040
40320
362880
3628800
39916800
479001600
6227020800
87178291200
1307674368000
20922789888000
355687428096000
6402373705728000
121645100408832000
2432902008176640000
-4249290049419214848
-1250660718674968576
8128291617894825984
-7835185981329244160
7034535277573963776
-1569523520172457984
-5483646897237262336
-5968160532966932480
-7055958792655077376
-8764578968847253504
4999213071378415616
-6045878379276664832

是否看到负数？那是您的（长整数）增长太多并从其最低可能值重新开始的时候。了解溢出；

顺便说一句，即使我将long替换为long long（并应用适当的printf格式%lld），我（在计算机上）也得到相同的结果。 >
计算此类内容非常复杂，并且需要大量有关数字如何存储在计算机中以及数值数学和数学分析的知识。我不认为自己是专家，所以这个问题可能比我写的要多。但是，我的解决方案似乎与C语言在其64位计算机上使用expl函数计算的结果一致，该函数使用gcc 4.7.2 20120921编译。

#2 楼

您的算法使用的中间数非常大和非常小，因此计算机很难使用它们。例如32！ ≈2.6313×1035，远超过LONG_MAX（可能小到231-1-≈2×109）。

基本上，一旦n足够大（也许13，也许21），factorial(n)返回一个“随机”数字。由于您的result是long而不是unsigned long，因此大多数溢出结果将是非常大的正数或非常大的负数，因此它们的倒数通常将平均为零，但并不能保证这一点。

所以，什么是好的算法？您需要避免使用很小或很大的中间值的东西。该问题被称为“堆栈溢出”问题3028312。该答案链接到描述许多计算算法的论文。本文的第9节有一个有效的算法（尽管它是代码跟踪的/神秘的）：


这是Xavier Gourdon编写的一个小型C程序，可以在您的计算机上计算9000的十进制数字电脑。对于π以及通过超几何级数定义的一些其他常数，存在相同类型的程序。

#include <stdio.h>
#define DIGITS 9000 /* decimal places (not including the '2') */
int main() {
    int N = DIGITS+9, a[DIGITS+9], x = 0;
    a[0] = 0;
    a[1] = 2;
    for (int n = 2; n < N; ++n) {
        a[n] = 1;
    }
    for ( ; N > 9; --N) {
        for (int n = N - 1; n > 0; --n) {
            a[n] = x % n;
            x = 10 * a[n-1] + x/n;
        }
        printf("%d", x);
    }
    return 0;
}

该程序（当有代码高尔夫时）具有117个字符。可以更改它以计算更多的数字（更改DIGITS的值）并且更快（将常数10更改为另一幂10并调整printf命令）。一个不太明显的问题是找到所使用的算法。

#3 楼

从技术上讲，问题不是下溢的问题，而是意义不大的问题。下溢是当数字太小而无法以当前格式存储时，将其替换为0.0。当将（相对）较大的数字添加到（相对）较小的数字时，就会失去重要性。回想一下，当添加分数时，必须将“十进制”（或二进制）点对齐。例如，在十进制中，我们可能会添加2.5x10 ^ 0和2.5x10 ^ -5，它们会扩展为将十进制对齐，如下所示...仅用4个精度执行此计算，答案就是2.500，而不是我们期望的2.500025。程序也会发生同样的事情，但是当然您有更多的位数可用，当然它以二进制而不是十进制出现，但是问题是相同的。