高斯积分和双除法

Question

为了好玩，我尝试使用级数展开来评估从 0 到 1 的高斯积分。出于这个原因，我写了一个阶乘函数，它在 20 以内都能很好地工作！（我检查过）然后我写了这个：

int main(){
    int n;
    long double result=0;
    for(n=0; n<=5; n++){
        if(n%2==0){
            result+=(((long double) 1/(long double)(factorial(n)*(2*n+1))));
        } else {
            result-=(((long double) 1/(long double)(factorial(n)*(2*n+1))));
        }
    }        
    printf("The Gaussian integral from 0 to 1 is %Lf\n", result);
}

这给了我一个奇怪的负数，显然甚至不接近。我怀疑问题出在演员表上，但我不知道是什么问题。有什么想法吗？这不是我尝试的第一件事。我尝试转换表达式中的任何内容并将显式转换放在开头，但它不起作用。

Answer 1

您正在使用 MinGW 编译器（Windows 的 gcc 端口），它与 long double 类型存在问题。这是由于 GCC 的 long double 实现与微软的 C 库之间的冲突。另见this question。

根据this question，定义__USE_MINGW_ANSI_STDIO 可以解决这个问题。如果没有，请改用double。

Answer 2

在(long double)(factorial(n)*(2*n+1) 中，乘法是整数乘法，如果factorial 的结果已经接近使用的整数类型的限制，第一个可能会溢出。

写入((long double)(factorial(n))*(2*n+1)，使第一个乘法是浮点乘法。

Answer 3

您几乎肯定会溢出您的整数类型。在 C 中，这在技术上是未定义的行为。

对于 32 位无符号整数，13！会溢出。在 64 位上，21！会溢出。

如果您使用浮点 double 类型或像 __uint128 之类的扩展名（我认为最多可以给您 34 个！），如果您的编译器支持，您的算法将存活更长时间。

您遇到的另一个问题是您正在逐步将减少大小的条款添加到您的总数中。在使用浮点类型时，这绝不是一个好主意。如果您以相反的顺序运行for 循环，那么结果会更准确。