/* main.c:	Main "Newton" program.
 *
 * Written by Daniel Barrett.  100% PUBLIC DOMAIN. */

#include "decl.h"
	
main(argc, argv)
int argc; char *argv[];
{
	complex polynomial[MAX_DEGREE+1];	/* Our polynomial.   */
	int degree;				/* Its degree.       */
	double epsilon;				/* Desired accuracy. */

#ifdef UNIX
	srand48((long)123);
#endif


	epsilon	= 0.0;
	degree	= 0;

	Version();
	InitPoly(polynomial, MAX_DEGREE+1);
	ReadPoly(polynomial, &degree, &epsilon);
	PrintPoly(polynomial, degree);
	NewtonMethod(polynomial, degree, epsilon);
}


NewtonMethod(poly, degree, epsilon)
/* Find the roots using Newton's Method. */
complex poly[];
int degree;
double epsilon;
{
	int numRoots=0, iterations;
	complex zero, oldGuess, newGuess, polyCopy[MAX_DEGREE+1];
	BOOL ErrorTooBig();

	AssignComplex(&zero, 0.0, 0.0);
	while (numRoots < degree-1) {		/* For each root...   */
		CopyPoly(poly, polyCopy, degree);
		Guess(&oldGuess);		/* Guess its value... */
		CopyComplex(oldGuess, &newGuess);
		iterations = MAX_ITERATIONS;
		do {				/* Refine the guess...*/
			CopyPoly(poly, polyCopy, degree);
			CopyComplex(newGuess, &oldGuess);
			SynthDivide(polyCopy, oldGuess, degree-numRoots);
			SynthDivide(polyCopy, oldGuess, (degree-1)-numRoots);
			Iterate(polyCopy, oldGuess, degree-numRoots, &newGuess);
		}while (ErrorTooBig(oldGuess, newGuess, epsilon)
		    &&  --iterations);		/* ...until convergence. */
		if (!iterations) {
			printf("Root didn't converge after %d iterations.\n",
				MAX_ITERATIONS);
			printf("Exiting!\n");
			exit(20);
		}
		SynthDivide(poly, newGuess, degree-numRoots);
		PrintRoot(newGuess, ++numRoots);
	}

	Divide(poly[1], poly[0], &poly[1]);	/* Get the last root. */
	Subtract(zero, poly[1], &newGuess);
	PrintRoot(newGuess, degree);
}