#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <strings.h>

#include "newton.h"

#define EPS 1e-6

int
close_to(struct z a, struct z b, double eps)
{
    double d1 = (a.re - b.re) / eps, d2 = (a.im - b.im) / eps;
	return (fabs(d1) < 1.0 && fabs(d2) < 1.0);
	return (d1 * d1 + d2 * d2) < 1.0;
}

void
newton(struct pngdata *pd)
{
	struct z z, dz, a = { -1.5, -1.0 }, b = { 1.5, 1.0 };
	double scale = 1.0;
	int i, j;
    struct roots roots;

    bzero((void *)&roots, sizeof(roots));

	a.re *= scale; a.im *= scale;
	b.re *= scale; b.im *= scale;

	if (pd->width <= 1 || pd->height <= 1) {
		fprintf(stderr, "Image dimensions are too small (%d x %d);\n"
		  "expecting at least 2 pixels.\n", pd->width, pd->height);
		exit (1);
	}

	// Working on rectangle (-a, -b) -- (a, b) at C
	dz.re = (b.re - a.re) / (pd->width - 1.0);
	dz.im = (b.im - a.im) / (pd->height - 1.0);

	for (i = 0, z.im = a.im; i < pd->height; i++, z.im += dz.im)
		for (j = 0, z.re = a.re; j < pd->width; j++, z.re += dz.re) {
			png_byte *pixel = pd->rows[i] + 3 * j;
			struct z p;
			int idx;
			unsigned long iters = newt_iter(z, &p, &roots, EPS);

			pixel[0] = pixel[1] = pixel[2] = 255;

			#if 1
			if (iters >= INFTY)
				continue;
			#endif

			for (idx = 0; idx < roots.n; idx++)
                if (close_to(p, roots.z[idx], EPS))
                    break;
            if (idx == roots.n || close_to(z, roots.z[idx], 1e-2)) {
				pixel[0] = pixel[1] = pixel[2] = 0;
				continue;
			}

			pixel[0] = pixel[1] = pixel[2] = 0;
			pixel[idx] = 255 * (1.0 - pow((double)iters / INFTY, 0.3));
		}
}