page@swan.ulowell.edu (Bob Page) (10/29/88)
Submitted-by: barrett@cs.jhu.edu (Dan Barrett)
Posting-number: Volume 2, Issue 37
Archive-name: applications/newton.1
Given a polynomial equation F(x) = 0, Newton estimates the roots of
the equation. The program uses an algorithm known as "Newton's
Method". You can read about the specifics in any college textbook on
Numerical Analysis.
# This is a shell archive. Remove anything before this line
# then unpack it by saving it in a file and typing "sh file"
# (Files unpacked will be owned by you and have default permissions).
# This archive contains the following files:
# ./Makefile
# ./Makefile.unix
# ./README
# ./complex.c
# ./decl.h
# ./ds.c
# ./io.c
# ./main.c
# ./strtok.c
#
if `test ! -s Makefile`
then
echo "writing Makefile"
cat > Makefile << '\Rogue\Monster\'
############################################################################
# Makefile for the "Newton" program, using Manx C V3.6a.
#
# Written by Daniel Barrett. 100% PUBLIC DOMAIN.
############################################################################
LFLAGS=-g +Q
SRC=main.c complex.c io.c ds.c strtok.c
OBJ=main.o complex.o io.o ds.o strtok.o
PROG=Newton
# Use ma32.lib for IEEE double precision floating point arithmetic.
CFLAGS=+L -n +Iinclude.pre +fi
LIBS=-lma32 -lc32
# Or, Use m32.lib for faster, less precise arithmetic.
#CFLAGS=+L -n +Iinclude.pre
#LIBS=-lm32 -lc32
# Or, use 68881 library for double precision arithmetic.
# Program crashed when I tried to run it after this compilation.
#CFLAGS=+L -n +Iinclude.pre +f8
#LIBS=-lm832 -lc32
amiga: include.pre $(PROG)
include.pre: decl.h
cc +L -n +Hinclude.pre decl.h
delete decl.o
$(PROG): $(OBJ)
ln $(LFLAGS) $(OBJ) -o $(PROG) $(LIBS)
$(OBJ): decl.h
clean:
delete \#?.o $(PROG) $(PROG).dbg include.pre
\Rogue\Monster\
else
echo "will not over write Makefile"
fi
if [ `wc -c Makefile | awk '{printf $1}'` -ne 994 ]
then
echo `wc -c Makefile | awk '{print "Got " $1 ", Expected " 994}'`
fi
if `test ! -s Makefile.unix`
then
echo "writing Makefile.unix"
cat > Makefile.unix << '\Rogue\Monster\'
############################################################################
# Makefile for the "Newton" program, UNIX 4.2BSD.
#
# Written by Daniel Barrett. 100% PUBLIC DOMAIN.
############################################################################
SRC=main.c complex.c io.c ds.c
OBJ=main.o complex.o io.o ds.o
LIBS=-lm
PROG=newton
CFLAGS=-DUNIX -g
all: $(PROG)
$(PROG): $(OBJ)
cc $(OBJ) -o $(PROG) $(LIBS)
$(OBJ): decl.h
clean:
rm -f *.o core
\Rogue\Monster\
else
echo "will not over write Makefile.unix"
fi
if [ `wc -c Makefile.unix | awk '{printf $1}'` -ne 466 ]
then
echo `wc -c Makefile.unix | awk '{print "Got " $1 ", Expected " 466}'`
fi
if `test ! -s README`
then
echo "writing README"
cat > README << '\Rogue\Monster\'
########################################################################
# Newton: Estimate the roots of a polynomial equation.
# Author: Dan Barrett, barrett@cs.jhu.edu (ARPAnet).
# 100% PUBLIC DOMAIN, both source code & executable.
#
# Written for the Commodore Amiga, September 1988.
# Runs from the CLI only.
# Also compiles & runs under Berkeley UNIX 4.2.
########################################################################
What is Newton?
---------------
Given a polynomial equation F(x) = 0, Newton estimates the
roots of the equation. The program uses an algorithm known as
"Newton's Method". You can read about the specifics in any college
textbook on Numerical Analysis.
How is it used?
---------------
From the CLI, type "newton", followed by <RETURN>.
You will be prompted for the degree of the polynomial equation.
The degree is the largest exponent in the equation. The maximum
allowable degree is 20. You can change this in the source code
by changing the value of MAX_DEGREE in decl.h.
Next, you are prompted for the "desired accuracy". Since
Newton's Method is a "numerical" method, the answers you get are
only estimates of the actual value.
The accuracy is a measure of how close an estimate must
be before it is considered "correct". In mathematical language,
the Euclidean distance (in the complex plane) between the current
estimate and the previous estimate must not exceed the "accuracy".
In the source code, the variable "epsilon" represents the
accuracy. See also the function ErrorTooBig(), which actually tests
the accuracy of the latest estimate.
To choose the default accuracy, simply press <RETURN>. Otherwise,
type in your desired accuracy. The smaller the value, the more accurate
your result (and the longer the program will run).
Finally, you are prompted for your polynomial coefficients.
For example, when you are asked:
X^5 coefficient?:
you should type in the coefficient of your X-to-the-5th-power term.
Coefficients may be real or complex numbers. If the coefficient
is real, simply type it in and press <RETURN>. If the coefficient is
complex, type in its real and imaginary parts, separated by spaces
and/or tabs, and press <RETURN>.
"Real" example:
X^5 coefficient?: 6.33 <RETURN>
"Complex" example, for the complex value -5 + 0.43i :
X^5 coefficient?: -5 0.43 <RETURN>
Note that you must type in EVERY coefficient, even if it is a zero.
While you are entering coefficients, you can type "H", "h", or "?"
for a brief "help" message.
Now What Happens?
-----------------
First, "Newton" will display your polynomial, as you typed
it in.
Then, "Newton" will calculate the roots of your polynomial
equation.
This calculation might take some time (a few minutes), depending
on the degree of the equation and the desired accuracy. Higher accuracy
takes more time, as you might guess. Sometimes, the calculation takes
only a few seconds.
Since "Newton" uses a numerical algorithm, there is the possibility
that it will not find a solution. If "Newton" cannot find the value of
a root after 10,000 iterations, it will report failure. (This number
of iterations is set in decl.h as MAX_ITERATIONS... feel free to change
it.)
So, if "Newton" seems to be sitting there not doing anything,
don't panic. After a few minutes, it will either find a solution
or quit automatically.
Compiling Information
---------------------
The enclosed "Makefile" is for use with Aztec Manx C, V3.6a.
"Newton" probably compiles & runs under Lattice C, though I can't
test it.
Math (floating point) calculations can be handled in three
different ways. Edit the Makefile and uncomment the desired
CFLAGS and LIBS parameters. They are explained in the Makefile itself.
See also the "New Options For CC" page (cc.ap.1, V3.4a) in the Manx
C manual.
The enclosed "Makefile.unix" is for compiling "Newton" under
UNIX (Berkeley UNIX 4.2 or Ultrix). It does not use strtok.c; I am
assuming that your UNIX has the strtok() function provided.
\Rogue\Monster\
else
echo "will not over write README"
fi
if [ `wc -c README | awk '{printf $1}'` -ne 4031 ]
then
echo `wc -c README | awk '{print "Got " $1 ", Expected " 4031}'`
fi
if `test ! -s complex.c`
then
echo "writing complex.c"
cat > complex.c << '\Rogue\Monster\'
/* complex.c: Complex arithmetic routines.
*
* Written by Daniel Barrett. 100% PUBLIC DOMAIN. */
#include "decl.h"
/* Complex arithmetic functions Add(), Subtract(), Multiply() and Divide()
* perform as their names indicate. They perform their operation on their
* first 2 arguments, and return the result in the third argument. */
Add(a, b, sum)
complex a, b, *sum;
{
sum->n[REAL] = a.n[REAL] + b.n[REAL];
sum->n[IMAG] = a.n[IMAG] + b.n[IMAG];
}
Subtract(a, b, diff)
complex a, b, *diff;
{
diff->n[REAL] = a.n[REAL] - b.n[REAL];
diff->n[IMAG] = a.n[IMAG] - b.n[IMAG];
}
Multiply(a, b, prod)
complex a, b, *prod;
{
prod->n[REAL] = (a.n[REAL] * b.n[REAL]) - (a.n[IMAG] * b.n[IMAG]);
prod->n[IMAG] = (a.n[REAL] * b.n[IMAG]) + (a.n[IMAG] * b.n[REAL]);
}
Divide(a, b, quot)
complex a, b, *quot;
{
double denom;
denom = SQR(b.n[REAL]) + SQR(b.n[IMAG]);
if (denom == 0.0)
printf("Divide by zero error!\n"), exit(20);
quot->n[REAL] = ((a.n[REAL] * b.n[REAL]) + (a.n[IMAG] * b.n[IMAG]))
/ denom;
quot->n[IMAG] = ((a.n[IMAG] * b.n[REAL]) - (a.n[REAL] * b.n[IMAG]))
/ denom;
}
SynthDivide(poly, comp, stop)
/* Computes the synthetic division of the polynomial poly[] by (z-comp),
* where "z" is assumed the complex variable in our polynomial. */
complex poly[], comp;
int stop;
{
int i;
complex prod;
for (i=1; i<=stop; i++) {
Multiply(poly[i-1], comp, &prod);
Add(poly[i], prod, &poly[i]);
}
}
Iterate(poly, zOld, n, zNew)
/* One iteration of Newton's method. "zOld" is the old guess of the value
* of a root of poly[]. "zNew" is the newly calculated guess. */
complex poly[], zOld;
int n;
complex *zNew;
{
complex quotient;
Divide(poly[n], poly[n-1], "ient);
Subtract(zOld, quotient, zNew);
}
Guess(z)
/* A random complex number, 50% chance of being negative. */
complex *z;
{
#ifdef UNIX
#define ran drand48
#endif
double ran();
z->n[REAL] = ran();
z->n[IMAG] = ran();
z->n[REAL] = (ran() < 0.50) ? z->n[REAL] : -(z->n[REAL]);
z->n[IMAG] = (ran() < 0.50) ? z->n[IMAG] : -(z->n[IMAG]);
}
BOOL ErrorTooBig(y, z, eps)
/* Compute the Euclidean distance between y and z in the complex plane.
* This is just our good old friend, the "distance formula". (Add the
* squares of y and z, and take the square root.) Is the distance less
* than epsilon? If so, the error is allowable; else, it's too big. */
complex y, z;
double eps;
{
complex diff;
double sqrt();
Subtract(y, z, &diff);
return( sqrt(SQR(diff.n[REAL]) + SQR(diff.n[IMAG])) < eps
? FALSE
: TRUE);
}
\Rogue\Monster\
else
echo "will not over write complex.c"
fi
if [ `wc -c complex.c | awk '{printf $1}'` -ne 2539 ]
then
echo `wc -c complex.c | awk '{print "Got " $1 ", Expected " 2539}'`
fi
if `test ! -s decl.h`
then
echo "writing decl.h"
cat > decl.h << '\Rogue\Monster\'
/* decl.h: Declarations for the "Newton" program.
*
* Written by Daniel Barrett. 100% PUBLIC DOMAIN. */
#include <stdio.h>
#include <math.h>
#ifdef UNIX
#define BOOL short
#define TRUE 1
#define FALSE 0
#else
#include <exec/types.h>
#endif
#define MAX_DEGREE 20 /* Maximum degree of equation. */
#define MAX_ITERATIONS 10000
#define EPSILON_DEFAULT 0.0000001
#define EQUAL !strcmp
#define SQR(x) ((x)*(x))
struct ComplexNumber { /* One complex number. */
double n[2]; /* Real & imaginary parts. */
};
typedef struct ComplexNumber complex;
#define REAL 0 /* 1st element of ComplexNumber. */
#define IMAG 1 /* 2nd element of ComplexNumber. */
/* If C is of type "complex", then:
*
* C.n[REAL] is the real part.
* C.n[IMAG] is the imaginary part.
*
*/
\Rogue\Monster\
else
echo "will not over write decl.h"
fi
if [ `wc -c decl.h | awk '{printf $1}'` -ne 782 ]
then
echo `wc -c decl.h | awk '{print "Got " $1 ", Expected " 782}'`
fi
if `test ! -s ds.c`
then
echo "writing ds.c"
cat > ds.c << '\Rogue\Monster\'
/* ds.c: Data structure routines.
*
* Written by Daniel Barrett. 100% PUBLIC DOMAIN. */
#include "decl.h"
InitPoly(poly, n)
/* Initialize the polynomial to all zeroes. */
complex poly[];
int n;
{
int i;
for (i=0; i<n; i++)
AssignComplex(&poly[i], 0.0, 0.0);
}
AssignComplex(comp, realPart, imagPart)
/* Assign the real & imaginary parts to a complex number. */
complex *comp;
double realPart, imagPart;
{
comp->n[REAL] = realPart;
comp->n[IMAG] = imagPart;
}
CopyPoly(poly1, poly2, degree)
/* Copy poly1 into poly2. */
complex poly1[], poly2[];
int degree;
{
int i;
for (i=0; i<=degree; i++)
CopyComplex(poly1[i], &poly2[i]);
}
CopyComplex(c1, c2)
/* Copy complex number c1 into c2. */
complex c1, *c2;
{
AssignComplex(c2, c1.n[REAL], c1.n[IMAG]);
}
\Rogue\Monster\
else
echo "will not over write ds.c"
fi
if [ `wc -c ds.c | awk '{printf $1}'` -ne 775 ]
then
echo `wc -c ds.c | awk '{print "Got " $1 ", Expected " 775}'`
fi
if `test ! -s io.c`
then
echo "writing io.c"
cat > io.c << '\Rogue\Monster\'
/* io.c: Input/Output functions.
*
* Written by Daniel Barrett. 100% PUBLIC DOMAIN. */
#include "decl.h"
PrintComplex(comp)
/* Print a complex number. */
complex comp;
{
double fabs();
if (comp.n[REAL] >= 0)
printf(" ");
printf("%f%c%fi", comp.n[REAL],
(comp.n[IMAG] < 0.0) ? '-' : '+',
fabs(comp.n[IMAG]));
}
PrintPoly(poly, n)
/* Print the polynomial in a nice format. */
complex poly[];
int n;
{
int i;
printf("\nYour polynomial is:\n");
printf("--------------------\n");
for (i=0; i<n-1; i++) {
PrintComplex(poly[i]);
printf(" * X^%d +\n", n-i);
}
PrintComplex(poly[n-1]);
printf(" * X +\n");
PrintComplex(poly[n]);
printf("\n\n");
}
ReadPoly(poly, degree, epsilon)
/* Read all data from the user: the degree of the polynomial, the
* accuracy desired, and the polynomial coefficients. */
complex poly[];
int *degree;
double *epsilon;
{
ReadDegree(degree);
ReadEpsilon(epsilon);
TellUserWhatToDo();
ReadCoefficients(poly, *degree);
}
ReadDegree(degree)
/* Prompt user for the degree of the polynomial. Read it. */
int *degree;
{
char buf[BUFSIZ];
*degree = 0;
do {
printf("Degree of polynomial? (1..%d): ", MAX_DEGREE);
gets(buf);
*degree = atoi(buf);
}while (*degree <= 0 || *degree > MAX_DEGREE);
}
ReadEpsilon(epsilon)
/* Prompt user for epsilon, the desired accuracy. Read it. If no
* data, set epsilon to the default value. */
double *epsilon;
{
char buf[BUFSIZ];
double atof();
do {
*epsilon = EPSILON_DEFAULT;
printf("Accuracy desired? [<RETURN> = %1.7f]: ", *epsilon);
gets(buf);
if (strlen(buf))
*epsilon = atof(buf);
}while (*epsilon <= 0.0);
}
ReadCoefficients(poly, degree)
/* Read in the coefficients of the polynomial "poly". NOTE that
* the array "poly" is backwards... poly[0] is the coefficient
* of X^degree, and so on.
* I use "strtok()", my version of the UNIX function, to get the input. */
complex poly[];
int degree;
{
char buf[BUFSIZ], *tokeReal=NULL, *tokeImag=NULL, *strtok();
static char separators[] = " \t\n\r";
register int i;
double atof();
for (i=0; i <= degree; i++) {
Again: printf("X^%d coefficient: ", degree-i);
gets(buf);
tokeReal = strtok(buf, separators);
if (!tokeReal)
goto Again;
else if (NeedsHelp(tokeReal)) {
Help();
goto Again;
}
tokeImag = strtok((char *)NULL, separators);
AssignComplex(&poly[i], atof(tokeReal),
tokeImag ? atof(tokeImag) : 0.0);
}
}
PrintRoot(root, i)
/* Print the value of root "i" of the polynomial. */
complex root;
int i;
{
printf("Root %d: ", i);
if (i<10)
printf(" ");
PrintComplex(root);
printf("\n");
}
TellUserWhatToDo()
{
printf("Enter the coefficients of your polynomial. ");
printf("Type \"H\" for help.\n");
}
static char *helpMessage[] = {
"",
"Enter the coefficients of your polynomial, one coefficient per line.",
"Note that \"X^n\" means \"X to the nth power\".",
"If the coefficient is REAL, simply enter it and press <RETURN>.",
"If the coefficient is COMPLEX, enter its real and imaginary terms,",
" separated by a space. (Do not use the letter \"i\".)",
"",
NULL
};
Help()
/* Print the above help message. */
{
char **s = helpMessage;
do
puts(*s);
while (*(s++));
}
NeedsHelp(s)
/* Does the user need help? */
char *s;
{
return(EQUAL(s, "H") | EQUAL(s, "h") | EQUAL(s, "?"));
}
Version()
{
printf("[33mNewton V1.0[0m");
printf(" by Daniel Barrett. 100%% PUBLIC DOMAIN.\n");
printf("Estimate the roots of a polynomial numerically.\n");
printf("Type ^C <RETURN> to abort this program.\n\n");
}
\Rogue\Monster\
else
echo "will not over write io.c"
fi
if [ `wc -c io.c | awk '{printf $1}'` -ne 3554 ]
then
echo `wc -c io.c | awk '{print "Got " $1 ", Expected " 3554}'`
fi
if `test ! -s main.c`
then
echo "writing main.c"
cat > main.c << '\Rogue\Monster\'
/* 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, °ree, &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);
}
\Rogue\Monster\
else
echo "will not over write main.c"
fi
if [ `wc -c main.c | awk '{printf $1}'` -ne 1760 ]
then
echo `wc -c main.c | awk '{print "Got " $1 ", Expected " 1760}'`
fi
if `test ! -s strtok.c`
then
echo "writing strtok.c"
cat > strtok.c << '\Rogue\Monster\'
/* strtok.c: 3 very useful string-handling functions.
* You don't need this file if you are using Lattice C.
* Manx C does not provide them.
*
* Written by Daniel Barrett. 100% PUBLIC DOMAIN. */
#include "decl.h"
#define STRING_END '\0'
/* Return the next string "token" in "buf". Tokens are substrings
* separated by any characters in "separators". */
char *strtok(buf, separators)
char *buf, *separators;
{
register char *token, *end; /* Start and end of token. */
extern char *strpbrk();
static char *fromLastTime;
if (token = buf ? buf : fromLastTime) {
token += strspn(token, separators); /* Find token! */
if (*token == STRING_END)
return(NULL);
fromLastTime = ((end = strpbrk(token,separators))
? &end[1]
: NULL);
*end = STRING_END; /* Cut it short! */
}
return(token);
}
/* Return a pointer to the first occurance in "str" of any character
* in "charset". */
char *strpbrk(str, charset)
register char *str, *charset;
{
extern char *index();
while (*str && (!index(charset, *str)))
str++;
return(*str ? str : NULL);
}
/* Return the number of characters from "charset" that are at the BEGINNING
* of string "str".
*/
int strspn(str, charset)
register char *str, *charset;
{
register char *s;
s = str;
while (index(charset, *s))
s++;
return(s - str);
}
\Rogue\Monster\
else
echo "will not over write strtok.c"
fi
if [ `wc -c strtok.c | awk '{printf $1}'` -ne 1322 ]
then
echo `wc -c strtok.c | awk '{print "Got " $1 ", Expected " 1322}'`
fi
echo "Finished archive 1 of 1"
# if you want to concatenate archives, remove anything after this line
exit
--
Bob Page, U of Lowell CS Dept. page@swan.ulowell.edu ulowell!page
Have five nice days.