Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!usc!wuarchive!udel!haven!adm!news
From: capehart@nevada.edu (Anne Racel)
Newsgroups: comp.lang.pascal
Subject: Re: Roots of polynomials
Message-ID: <25249@adm.brl.mil>
Date: 12 Dec 90 02:40:36 GMT
Sender: news@adm.brl.mil
Lines: 95


=> 	Finding roots for a quadratic was easy enough.  However,
=> 	if you could give me pointers to finding roots to higher
=> 	degree polynomials, I'd be interested in hearing from you.
=> 	(no, this is not a lab assignment... :)

   What I used, courtesy of my chem instructor (was transferring a program
for pH from Fortran to Pascal, was the Newton-Reismann (sp?) formula.
from calculus.  Must admit, there is some bug that hopefully some guru out
there can figure out, that causes the thing to take 20 mintues for a
wrongly rounded anser.  But here, essentially, is it:
 
 
function  NEWTON(TERMS   : POLY  ;
                 OLD     : real): real;
 
var
   ORDER : integer   ;                   { order of polynomial expression }
   H,I     : integer ;                   { counter                        }
   DERIV   : real    ;                   { derivative                     }
   Y       : real    ;                   { haven't a clue                 }
   CRIT    : real    ;                   { criteria for exiting loop      }
   DUM     : real    ;                   { haven't a clue                 }
   NEW     : real    ;                   { new value for H+ concentration }
 
begin
 
   if (TERMS[1]) <> 0 then
      begin
         Y     := 0;
         ORDER := 12;
         DERIV := 0 ;
         repeat                               { find order of the polynomial }
            ORDER := ORDER - 1
         until TERMS[ORDER] > 0;
         { end repeat }
         if ORDER > 1 then
            repeat
               for I := 1 to ORDER do          {evaluate polynomial}
                  begin
                     Y:= Y + TERMS[I] * EXP(I * LN(OLD));
                     DERIV := DERIV + ((I-1) * TERMS[I] * EXP((I-2)*LN(OLD)));
                  end;
               NEW := OLD - (Y/DERIV);
               CRIT  := ABS(1.0e6 * ((OLD-NEW)/OLD));       {check whether old value and new value close }
               OLD := NEW;
               Y := 0;
               DERIV := 0
             until CRIT <= 1
         else
            NEW := OLD
       end;
       NEWTON := NEW
end;
 

..... and for those familiar with Fortran, here's that version, which DOES
work:

$STORAGE:2
      SUBROUTINE NEWTON (POLY,X,NEWX)
      INTEGER H, I
      INTEGER*4 CRIT
      REAL*8 POLY, DERIV, X, Y, NEWX, DUM
      DIMENSION POLY(11)
C     FIND THE ORDER OF THE POLYNOMIAL EXPRESSION STORED IN THE MATRIX
      DO 5, H=11, 1, -1
      IF (POLY(H)) 10, 5, 10
5     CONTINUE
10    DO 20, I=1, H
      Y = Y + POLY(I)*X**(I-1)
      IF (I.LE.1) GOTO 20
      DERIV = DERIV + (I-1)*POLY(I)*X**(I-2)
20    CONTINUE
      NEWX = X - (Y/DERIV)
      DUM = 1000000.*((X-NEWX)/X)
      CRIT = ABS(INT(DUM))
      IF (CRIT.LE.1) GOTO 30
      X = NEWX
      Y = 0
      DERIV = 0
      GOTO 10
30    RETURN
      END

 
   ============================================================= 
    Anne Racel                Internet:  capehart@nevada.edu
    University of Nevada,   Compuserve:  72105,1105
    Las Vegas                   Bitnet:  capehart@unsvax.bitnet

    "I'm very certain, Oz, that you gave me the best brains in the world
    for I can think with them day and night, when all over brains are
    fast asleep."
			Scarecrow in "Dorothy and the Wizard in Oz"