复数 (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 的限制，后者将 return 1，这对于任何限制都是正确的，而前者将 return 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

输入的差异 −2 + 0 i

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

z = -2 + 0 i.
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 = z*z + c替换z = cpowl(z,2)+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

Answer 2

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

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

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

zi 的计算是使用 zr 的 new 值，而不是当前值。您需要将这两个计算的结果保存在临时变量中，然后将它们分配回 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 来避免。

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

Complex numbers (complex.h) and apparent lag of precision

c

complex-numbers

输入的差异 −2 + 0 i

输入的差异 −1.3 + 0.1 i