c - 为什么C程序第10次循环结束时循环语句会出错？

Question

我正在尝试计算序列的近似值，并将精度设置为 1E-6。

int i=1;
    float x;
    scanf("%f",&x);
    double sx = x,temp;
    do{
        temp = (pow(x,((2*i)+1)))/((2*i+1)*(fact(i))); //fact() is a function counting factorial
        if((i+2)%2==1) sx-=temp;                        
        else sx+=temp;
        printf("%lf %lf\n",temp,sx); 
        i++;
    }while(temp>1E-6);
    printf("%lf",sx);

当输入小于2时，我得到预期的输出。但是当输入达到3时，出现负数。是因为堆栈溢出吗？但事实(10)并没有那么大，不是吗？我怎样才能得到正确的结果？

Answer 1

你问的是不是10!太大了，但是如果您在值溢出时观察 i，您会看到 i == 12。 temp then 的表达式为:

3²⁵/25 x 12!

子表达式 25 x 12!值为 11975040000，如果它是一个整数表达式，则需要 34 位，如果 fact() 是一个整数函数，则需要 34 位。如果是32位整数函数就会溢出。

如果将 fact() 更改为无符号 64 位类型，它会更进一步，但这仅适用于 19!并且除数子表达式在收敛之前在 i=18 (29 x 18!) 处溢出。

使用以下阶乘函数实现，它将收敛于 i = 3:

double fact(int n)
{
    if (n >= 1)
        return n * fact(n - 1);
    else
        return 1;
}

结果为 0.886207。对于任意大的 x 值，您将需要任意精度或“bignum”数学库，而不是内置类型。

受到 Lee Crocker 评论的启发，我尝试了各种表达方式(在 GDB Online 进行了测试)

1) 基于我的双重事实( int ) 并减去你的无关括号；

double e = (2*i)+1 ;
temp = pow(x,e) / (e * fact(i));

x = 1 到 9 的结果:

x: 1.0  s(x): 0.746824                                                                                                                          
x: 2.0  s(x): 0.882081                                                                                                                          
x: 3.0  s(x): 0.886207                                                                                                                          
x: 4.0  s(x): 0.886227                                                                                                                          
x: 5.0  s(x): 0.886227                                                                                                                          
x: 6.0  s(x): 0.889677                                                                                                                          
x: 7.0  s(x): 2453.707036                                                                                                                       
x: 8.0  s(x): 1055556860.115990                                                                                                                 
x: 9.0  s(x): -nan 

2) 使用tgamma():

temp = pow(x,e) / (e * tgamma(i+1)) ; 

结果:

x: 1.0  s(x): 0.746824                                                                                                                          
x: 2.0  s(x): 0.882081                                                                                                                          
x: 3.0  s(x): 0.886207                                                                                                                          
x: 4.0  s(x): 0.886227                                                                                                                          
x: 5.0  s(x): 0.886226                                                                                                                          
x: 6.0  s(x): 0.882858                                                                                                                          
x: 7.0  s(x): -83.077451                                                                                                                        
x: 8.0  s(x): 3729162898.780624  

3) 使用lgamma():

temp = pow(x,e) / exp(log(e) + lgamma(i+1)) ;

结果:

x: 1.0  s(x): 0.746824                                                                                                                          
x: 2.0  s(x): 0.882081                                                                                                                          
x: 3.0  s(x): 0.886207                                                                                                                          
x: 4.0  s(x): 0.886227                                                                                                                          
x: 5.0  s(x): 0.886236                                                                                                                          
x: 6.0  s(x): 1.904565                                                                                                                          
x: 7.0  s(x): -322346.422518                                                                                                                    
x: 8.0  s(x): 5416428382.734000                                                                                                                 
x: 9.0  s(x): -nan       

对于x = 6.0，结果有所不同，但所有版本都在此处崩溃。我不清楚哪个更好 - 我怀疑 fact() 。标准库函数没有提供任何重大改进。

4)添加Lee的表情进行测试(如果我正确理解了他的答案):

temp = exp(e * log(x) - log(e) - lgamma(i+1)) ;

结果:

x: 1.0  s(x): 0.746824
x: 2.0  s(x): 0.882081
x: 3.0  s(x): 0.886207
x: 4.0  s(x): 0.886227
x: 5.0  s(x): 0.886245
x: 6.0  s(x): 2.892730
x: 7.0  s(x): 96596.597282
x: 8.0  s(x): -720878097610.696167
x: 9.0  s(x): -18109761159507625984.000000

仅供引用，我的测试代码:

double s( double x )
{
    double sx = x ;
    double temp = 1.0 ;

    for( int i = 1; temp > 1e-6; i++ ) 
    {
        double e = (2*i)+1 ;

        // Uncomment one of the following:
        //temp = pow(x,e) / (e * fact(i));
        //temp = pow(x,e) / (e * tgamma(i+1)) ; 
        //temp = pow(x,e) / exp(log(e) + lgamma(i+1)) ;
        //temp = exp(e * log(x) - log(e) - lgamma(i+1)) ;

        sx += (i % 2) == 0 ? temp : -temp ;  
    }

    return sx ;
}

int main()
{
    for( double x = 1.0; x < 10.0; x += 1.0 )
    {
        printf( "x: %.1f  s(x): %lf\n", x, s(x) ) ;
    }

    return 0;
}