Бассейны Ньютона для уравнения z^4 - 1 = 0

/*

    Бассейны Ньютона для уравнения z^4 - 1 = 0

*/

#include <math.h>
#include <stdio.h>
#include <string.h>



/* всякие функции, эмулирующие работу с комплексными числами */
struct Complex
    {
        float x;
        float y;
    };

typedef struct Complex   ComplexType;
typedef struct Complex * ComplexRef;

ComplexType Complex_number (float x, float y)
{
    ComplexType result;
    result.x = x;
    result.y = y;
    return result;
}

ComplexType Complex_sum (ComplexType u, ComplexType v)
{
    ComplexType result;
    result.x = u.x + v.x;
    result.y = u.y + v.y;
    return result;
}

ComplexType Complex_sub (ComplexType u, ComplexType v)
{
    ComplexType result;
    result.x = u.x - v.x;
    result.y = u.y - v.y;
    return result;
}

ComplexType Complex_prod (ComplexType u, ComplexType v)
{
    ComplexType result;
    result.x = u.x*v.x - u.y*v.y,
    result.y = u.x*v.y + u.y*v.x;
    return result;
}

float Complex_abs2 (ComplexType u)
{
    return u.x*u.x + u.y*u.y;
}

ComplexType Complex_conj (ComplexType u)
{
    return Complex_number (u.x, -u.y);
}

ComplexType Complex_div (ComplexType u, ComplexType v, int *error)
{
    *error = 0;
    float
        r = Complex_abs2(v);
    if (r == 0)
        *error = 1;
    ComplexType result = Complex_prod(u, Complex_conj(v));
    if (!(*error))
    {
        result.x /= r;
        result.y /= r;
    }
    return result;
}

ComplexType Complex_pown(ComplexType z, int n)
{
    ComplexType result = Complex_number(1, 0);
    while(n > 0)
    {
        if (n & 1)
            result = Complex_prod(result, z);
        n >>= 1;
        z = Complex_prod(z, z);
    }
    return result;
}
/* end of всякие функции, эмулирующие работу с комплексными числами */



/* функция, определяющая, к какому из четырех корней {1, -1, i, -i} уравнения z^4 - 1 = 0
   сходится итерационный процесс метода Ньютона от начальной точки z0 = x + iy */
int newton (float x, float y)
{
    ComplexType
        z = Complex_number(x, y),
        prev_z = Complex_number(1000, 1000);

    int error;

    /* пока очередное приближение не стало слишком большим (процесс расходится) и
    очередное приближение заметно отличается от предыдущего (процесс сошелся), считаем итерации
    */
    while ((Complex_abs2(z) < 1e6) && (Complex_abs2(Complex_sub(z, prev_z)) > 1e-6))
    {
        // запоминаем предыдущее приближение
        prev_z = z;

        // считаем новое приближение по формуле z - (z^4 - 1) / (4 z^3)
        z = Complex_sub(z, Complex_div(Complex_sub(Complex_pown(z, 4), Complex_number(1,0)), Complex_prod(Complex_number(4,0), Complex_pown(z, 3)), &error));

        // если на (4 z^3) поделить нельза, заплакали и вышли
        if (error)
            return 0;
    }
    /* end of пока очередное приближение ... считаем итерации*/

    /* анализируем результат: */

    // процесс разошелся
    if (Complex_abs2(z) > 1e6)
        return 0;

    // нас прибило к корню 1
    if (Complex_abs2(Complex_sub(z, Complex_number(1,0))) < 1e-3)
        return 1;

    // нас прибило к корню -1
    if (Complex_abs2(Complex_sub(z, Complex_number(-1,0))) < 1e-3)
        return 2;

    // нас прибило к корню i
    if (Complex_abs2(Complex_sub(z, Complex_number(0,1))) < 1e-3)
        return 3;

    // нас прибило к корню -i
    if (Complex_abs2(Complex_sub(z, Complex_number(0,-1))) < 1e-3)
        return 4;

    // нас занесло неведомо куда
    return 0;
}


int main(void)
{
    // задаем размеры картинки
    int Xpixels = 500, Ypixels = 500;

    // задаем размеры области комплексной плоскости
    float Xmin = -10.0, Xmax = 10.0, Ymin = -10.0, Ymax = 10.0;

    /* рисуем картинку */
    FILE *f = fopen("result.ppm", "w");
    fprintf(f, "P3\n");
    fprintf(f, "%d %d 255\n", Xpixels, Ypixels);

    for (int y = 0; y < Ypixels; y++)
    {
        for (int x = 0; x < Xpixels; x++)
        {
            switch (newton(Xmin + (Xmax - Xmin) / Xpixels * x, Ymin + (Ymax - Ymin) / Ypixels * y))
            {
                case 1: fprintf(f, "%d %d %d  ", 255,   0,   0); break;
                case 2: fprintf(f, "%d %d %d  ",   0, 255,   0); break;
                case 3: fprintf(f, "%d %d %d  ",   0,   0, 255); break;
                case 4: fprintf(f, "%d %d %d  ", 255, 255,   0); break;
                default:
                        fprintf(f, "%d %d %d  ",   0,   0,   0); break;
            }
        }
        fprintf(f, "\n");
    }

    fclose(f);
    /* end of рисуем картинку */

    return 0;
}