c - 为什么 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 都以给定的精度舍入到相同的数字。使用 long double 和上述函数，调试器显示前 2 次迭代:

    a = 1.5555555555555556
    a = 1.6444444444444444

虽然它应该是:

    a = 1.5555555555555556
    a = 1.653104925053533

此后程序无法更新值。我使用的线性插值方程是给定的数学方程的重新排列版本 here我使用的代码是我编写的 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() (来自 <stdlib.h> )具有签名 int abs(int); — 你从你的计算中得到整数。

你应该使用 long double fabsl(long double);来自 <math.h> .

您还应该使用 powl()而不是 pow() ( long double 与 double )和 roundl()而不是 roundf() ( long double 对比 float)。

换句话说，确保您使用的是正确的类型。

---

当您解决了类型问题后，您仍然会遇到收敛问题。如果您制作了 MCVE ( Minimal, Complete, Verifiable Example )，它会有所帮助。但是，这是我可以从您的问题中推断出的 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在给定的小数位数之内——或者类似的东西。