复数 (complex.h) 和明显的精度滞后

Question

我决定尝试一下 complex.h，但遇到了一个我认为非常奇怪的问题。

int mandelbrot(long double complex c, int lim)
{
    long double complex z = c;
    for(int i = 0; i < lim; ++i, z = cpowl(z,2)+c)
    {
        if(creall(z)*creall(z)+cimagl(z)*cimagl(z) > 4.0)
            return 0;
    }
    return 1;
}

int mandelbrot2(long double cr, long double ci, int lim)
{
    long double zr = cr;
    long double zi = ci;
    for(int i = 0; i < lim; ++i, zr = zr*zr-zi*zi+cr, zi = 2*zr*zi+ci)
    {
        if(zr*zr+zi*zi > 4.0)
            return 0;
    }
    return 1;
}

这些函数的行为不同。如果我们输入 -2.0+0.0i 和高于 17 的限制，后者将返回 1，这对于任何限制都是正确的，而前者将返回 0，至少在我的系统上是这样。 GCC 9.1.0，锐龙 2700x。

我这辈子都想不通这是怎么发生的。我的意思是，虽然我可能不完全理解 complex.h 在幕后是如何工作的，但对于这个特定的示例，结果应该像这样偏离是没有意义的。

---

在写作时，我注意到 cpowl(z,2)+c，并尝试将其更改为 z*z+c，这有所帮助，但是经过快速测试，我发现行为仍然不同。前任。 -1.3+0.1*I，lim=18。

我很想知道这是否特定于我的系统以及可能的原因是什么，尽管我完全清楚最相似的情况是我犯了一个错误，但是很遗憾，我找不到它.

--- 编辑---

最后是完整的代码，包括更改和修复。这两个函数现在似乎产生了相同的结果。

#include <stdio.h>
#include <complex.h>

int mandelbrot(long double complex c, int lim)
{
    long double complex z = c;
    for(int i = 0; i < lim; ++i, z = z*z+c)
    {
        if(creall(z)*creall(z)+cimagl(z)*cimagl(z) > 4.0)
            return 0;
    }
    return 1;
}

int mandelbrot2(long double cr, long double ci, int lim)
{
    long double zr = cr;
    long double zi = ci;
    long double tmp;
    for(int i = 0; i < lim; ++i)
    {
        if(zr*zr+zi*zi > 4.0) return 0;
        tmp = zi;
        zi = 2*zr*zi+ci;
        zr = zr*zr-tmp*tmp+cr;
    }
    return 1;
}

int main()
{
    long double complex c = -2.0+0.0*I;
    printf("%i\n",mandelbrot(c,100));
    printf("%i\n",mandelbrot2(-2.0,0.0,100));
    return 0;
}

cpowl() 仍然搞砸了，但我想如果我愿意，我可以创建自己的实现。

Answer 1

第二个函数是不正确的，而不是第一个。

在for的第三个子句的表达式中：

zr = zr*zr-zi*zi+cr, zi = 2*zr*zi+ci

zi 的计算是使用zr 的新 值，而不是当前值。您需要将这两个计算的结果保存在临时变量中，然后将它们分配回zr 和zi：

int mandelbrot2(long double cr, long double ci, int lim)
{
    long double zr = cr;
    long double zi = ci;
    for(int i = 0; i < lim; ++i)
    {   
        printf("i=%d, z=%Lf%+Lfi\n", i, zr, zi);
        if(zr*zr+zi*zi > 4.0)
            return 0;
        long double new_zr = zr*zr-zi*zi+cr;
        long double new_zi = 2*zr*zi+ci;
        zr = new_zr; 
        zi = new_zi;
    }
    return 1;
}

此外，使用cpowl 进行简单平方会导致不准确，在这种情况下可以通过简单地使用z*z 来避免。

Answer 2

输入差异 -2 + 0 i

cpowl 不准确。求幂是一个实现起来很复杂的函数，在其计算中可能会出现各种错误。在 macOS 10.14.6 上，mandelbrot 例程中的 z 在连续迭代中采用这些值：

z = -2 + 0 我。 z = 2 + 4.33681e-19 i。 z = 2 + 1.73472e-18 i。 z = 2 + 6.93889e-18 i。 z = 2 + 2.77556e-17 i。 z = 2 + 1.11022e-16 i。 z = 2 + 4.44089e-16 i。 z = 2 + 1.77636e-15 i。 z = 2 + 7.10543e-15 i。 z = 2 + 2.84217e-14 i。 z = 2 + 1.13687e-13 i。 z = 2 + 4.54747e-13 i。 z = 2 + 1.81899e-12 i。 z = 2 + 7.27596e-12 i。 z = 2 + 2.91038e-11 i。 z = 2 + 1.16415e-10 i。 z = 2 + 4.65661e-10 i。

因此，一旦出现初始错误，产生 2 + 4.33681•10⁻¹⁹ i，z 继续增长（正确地，由于数学，不仅仅是浮点错误) 直到它大到足以通过将其绝对值的平方与 4 进行比较的测试。（测试不会立即捕获超出部分，因为虚部的平方非常小，当添加到的平方时会丢失四舍五入真实的部分。）

相反，如果我们将z = cpowl(z,2)+c 替换为z = z*z + c，z 仍为 2（即 2 + 0i）。一般来说，z*z 中的操作也会遇到一些舍入错误，但没有cpowl 中的那么严重。

输入差异 -1.3 + 0.1 i

对于这个输入，差异是由于for循环的更新步骤中计算不正确造成的：

++i, zr = zr*zr-zi*zi+cr, zi = 2*zr*zi+ci

在计算zi 时使用zr 的新值。可以通过插入long double t;并将更新步骤更改为

来修复

++i, t = zr*zr - zi*zi + cr, zi = 2*zr*zi + ci, zr = t