#include <math.h>               /* sqrt() */
#include <complex>
#include "NewtonComputation.h"
using namespace std;

#define LEEWAY 0.001
#define ITERATIONS 32
#define CLOSE_ENOUGH(z, x, y) \
    (fabs(real(z) - (x)) < LEEWAY && fabs(imag(z) - (y)) < LEEWAY)

/*
 * This Computation class generates a design in the complex plane
 *  which shows which complex points go to which of the six roots
 *  of the equation z^6-1=0 using Newton's method of finding roots
 */
unsigned int NewtonComputation::compute(double x, double y)
{
    complex<double> rootGuess(x, y);
    static double halfRoot3 = sqrt(3) / 2.0;

    for (int iter = 0; iter < ITERATIONS; iter++)
    {
        /* percentage of color to illuminate based on current iteration */
        double pIlum = (double) (ITERATIONS - iter) / (double) ITERATIONS;
        /*
         * These if statements check to see if the complex number (the root guess)
         *  is within LEEWAY distance of a real root. If so, a unique color is returned
         *  reflecting which root it is close to and how many iterations it took
         *  for the root guess to get to that root
         */
        if (CLOSE_ENOUGH(rootGuess, 1, 0))
        {
            return ((unsigned int)(0xFF * pIlum)) << 16;
        }
        if (CLOSE_ENOUGH(rootGuess, -1, 0))
        {
            return (((unsigned int)(0xFF * pIlum)) << 16) + (((unsigned int)(0xFF * pIlum)) << 8);
        }
        if (CLOSE_ENOUGH(rootGuess, 0.5, halfRoot3))
        {
            return ((unsigned int)(0xFF * pIlum)) << 8;
        }
        if (CLOSE_ENOUGH(rootGuess, -0.5, halfRoot3))
        {
            return (unsigned int)(0x88 * pIlum);
        }
        if (CLOSE_ENOUGH(rootGuess, 0.5, -halfRoot3))
        {
            return (unsigned int)(0xFF * pIlum);
        }
        if (CLOSE_ENOUGH(rootGuess, -0.5, -halfRoot3))
        {
            return (((unsigned int)(0xFF * pIlum)) << 16) + (unsigned int)(0xFF * pIlum);
        }

        /* This expression evaluates the next complex number to be tested
         * for being close to a root. It uses Newton's method for finding
         * roots of equations according to the following recursive equation:
         *     x_n+1 = x_n - y_n / dy_n
         * --> x_n+1 = x_n - x_n^6 / 6x_n^5
         * More information can be found at:
         * http://www.willamette.edu/~sekino/fractal/fractal.htm
         */
        rootGuess -= (pow(rootGuess, 6.0) - 1.0) / (6.0 * pow(rootGuess, 5.0));
    }
    return 0;
}