67 lines
2.5 KiB
C++
67 lines
2.5 KiB
C++
|
|
#include <math.h> /* sqrt() */
|
|
#include <complex.h> /* complex, I, creal(), cimag() */
|
|
#include "NewtonComputation.h"
|
|
|
|
#define LEEWAY 0.001
|
|
#define ITERATIONS 32
|
|
#define CLOSE_ENOUGH(z, x, y) \
|
|
(fabs(creal(z) - (x)) < LEEWAY && fabs(cimag(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)
|
|
{
|
|
double complex rootGuess = x + I * 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 -= (cpow(rootGuess, 6.0) - 1.0) / (6.0 * cpow(rootGuess, 5.0));
|
|
}
|
|
}
|