复数 (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 来避免。