I decided to play around a bit with complex.h, and ran into what I consider a very curious problem.
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;
}
These functions do not behave the same. If we input -2.0+0.0i and a limit higher than 17, the latter will return 1, which is correct for any limit, while the former will return 0, at least on my system. GCC 9.1.0, Ryzen 2700x.
I cannot for the life of me figure out how this can happen. I mean while I may not entirely understand how complex.h works behind the scenes, for this particular example it makes no sense that the results should deviate like this.
While writing I notices the cpowl(z,2)+c, and tried to change it to z*z+c, which helped, however after a quick test, I found that the behavior still differ. Ex. -1.3+0.1*I, lim=18.
I'm curious to know if this is specific to my system and what the cause might be, though I'm perfectly aware that the most like scenario is me having made a mistake, but alas, I can't find it.
--- edit---
Finally, the complete code, including alterations and fixes. The two functions now seem to yield the same result.
#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() still messes things up, but I suppose if I wanted to, I could just create my own implementation.
cpowl is inaccurate. Exponentiation is a complicated function to implement, and a variety of errors likely arise in its computation. On macOS 10.14.6, z in the mandelbrot routine takes on these values in successive iterations:
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.
Thus, once the initial error is made, producing 2 + 4.33681•10−19 i, z continues to grow (correctly, as a result of mathematics, not just floating-point errors) until it is large enough to pass the test comparing the square of its absolute value to 4. (The test does not immediately capture the excess because the square of the imaginary part is so small it is lost in rounding when added to the square of the real part.)
In contrast, if we replace z = cpowl(z,2)+c with z = z*z + c, z remains 2 (that is, 2 + 0i). In general, the operations in z*z experience some rounding errors too, but not as badly as with cpowl.
For this input, the difference is caused by the incorrect calculation in the update step of the for loop:
++i, zr = zr*zr-zi*zi+cr, zi = 2*zr*zi+ci
That uses the new value of zr when calculating zi. It can be fixed by inserting long double t; and changing the update step to
++i, t = zr*zr - zi*zi + cr, zi = 2*zr*zi + ci, zr = t