为什么C为此算法赋予Python不同的值？

Question

我在C中写了一个简短的程序来执行线性插值,它迭代给出一个函数的根到给定的小数点数.到目前为止,对于函数f(x):

    long double f(long double x) {
        return (pow(x, 3) + 2 * x - 8);
    }

该程序无法收敛到1dp值.程序更新变量a和b,f(x)的根位于变量a和b之间,直到a和b都以给定的精度舍入到相同的数字.使用long double和上面的函数,调试器显示前2次迭代:

    a = 1.5555555555555556
    a = 1.6444444444444444

虽然应该是:

    a = 1.5555555555555556
    a = 1.653104925053533

在此之后,程序无法更新值.我正在使用的线性插值方程是这里给出的数学方法的重新排列版本,我使用的代码是我编写的python程序的C版本.尽管算法相同,为什么C实现会得到不同的值,我该如何解决？

好吧,我仍然掌握这一点,但希望下面我有一个最小的,完整的,可验证的例子:

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

long double a; long double b; long double c; // The values for interpolation
long double fofa; long double fofb; long double fofc; // The values f(a), f(b) and f(c)
const int dp = 1; // The number of decimal places to be accurate to

long double f(long double x) {
    return (pow(x, 3) + 2 * x - 8);
}

int main(void) {
    a = 1; b = 2;
    while(roundf(a * pow(10, dp))/pow(10, dp) != roundf(b * pow(10, dp))/pow(10, dp)) { // While a and b don't round to the same number...

        fofa = f(a); fofb = f(b); // Resolve the functions
        printf("So f(a) = %g, f(b) = %g\n", (double)fofa, (double)fofb); // Print the results

        c = (b * abs(fofa) + a * abs(fofb)) / (abs(fofb) + abs(fofa)); // Linear Interpolation

        fofc = f(c);

        if(fofc < 0) {
            a = c;
        }
        else if(fofc == 0) {
            a = c;
            break;
        }
        else {
            b = c;
        }

    }

    printf("The root is %g, to %d decimal places, after %d iterations.\n", (double)a, dp, i);
}

Answer 1

函数abs()(from <stdlib.h>)具有签名int abs(int);- 您将从计算中获得整数.

你应该用long double fabsl(long double);的<math.h>.

你也应该使用powl()而不是pow()(long doublevs double),而roundl()不是roundf()(long doublevs float).

换句话说,确保使用正确的类型.

当您修复类型问题时,仍然存在收敛问题.如果您生成了MCVE(最小,完整,可验证的示例),那将会有所帮助.但是,这是我可以从您的问题中推断出的MCVE:

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

static inline long double f(long double x)
{
    return(powl(x, 3) + 2 * x - 8);
}

int main(void)
{
    long double a = 1.0L;
    long double b = 2.0L;
    int dp = 6;

    while (roundl(a * powl(10, dp)) / powl(10, dp) != roundl(b * powl(10, dp)) / powl(10, dp))
    {
        long double fofa = f(a);
        long double fofb = f(b);
        long double c = (b * fabsl(fofa) + a * fabsl(fofb)) / (fabsl(fofb) + fabsl(fofa));
        long double fofc = f(c);
        printf("a = %+.10Lf, f(a) = %+.10Lf\n", a, fofa);
        printf("b = %+.10Lf, f(b) = %+.10Lf\n", b, fofb);
        printf("c = %+.10Lf, f(c) = %+.10Lf\n", c, fofc);
        putchar('\n');

        if (fofc < 0.0L)
        {
            a = c;
        }
        else if (fofc == 0.0L)
        {
            a = c;
            break;
        }
        else
        {
            b = c;
        }
    }
    printf("Result: a = %Lg\n", a);
    return 0;
}

我得到的输出是:

a = +1.0000000000, f(a) = -5.0000000000
b = +2.0000000000, f(b) = +4.0000000000
c = +1.5555555556, f(c) = -1.1248285322

a = +1.5555555556, f(a) = -1.1248285322
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6531049251, f(c) = -0.1762579238

a = +1.6531049251, f(a) = -0.1762579238
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6677455452, f(c) = -0.0258828049

a = +1.6677455452, f(a) = -0.0258828049
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6698816424, f(c) = -0.0037639074

a = +1.6698816424, f(a) = -0.0037639074
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6701919841, f(c) = -0.0005465735

a = +1.6701919841, f(a) = -0.0005465735
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702370440, f(c) = -0.0000793539

a = +1.6702370440, f(a) = -0.0000793539
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702435859, f(c) = -0.0000115206

a = +1.6702435859, f(a) = -0.0000115206
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702445357, f(c) = -0.0000016726

a = +1.6702445357, f(a) = -0.0000016726
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446735, f(c) = -0.0000002428

a = +1.6702446735, f(a) = -0.0000002428
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446936, f(c) = -0.0000000353

a = +1.6702446936, f(a) = -0.0000000353
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446965, f(c) = -0.0000000051

a = +1.6702446965, f(a) = -0.0000000051
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446969, f(c) = -0.0000000007

a = +1.6702446969, f(a) = -0.0000000007
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446970, f(c) = -0.0000000001

a = +1.6702446970, f(a) = -0.0000000001
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446970, f(c) = -0.0000000000

a = +1.6702446970, f(a) = -0.0000000000
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446970, f(c) = -0.0000000000

a = +1.6702446970, f(a) = -0.0000000000
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446970, f(c) = -0.0000000000

a = +1.6702446970, f(a) = -0.0000000000
b = +2.0000000000, f(b) = +4.0000000000
c = +1.6702446970, f(c) = -0.0000000000

无限循环的原因很明显; 之间的差a和b不小的部分.您需要查看收敛标准.它可能应该fofc与0.0给定的小数位数之间进行比较- 或者沿着这些行进行比较.