CTITLEMPOLRT -- COMPUTES THE REAL AND COMPLEX ROOTS OF POLYNOMIAL 00000100 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR IBM 00000200 CA DESIGNER IBM 00000300 CA LANGUAGE FORTRAN H 00000400 CA SYSTEM S/370 00000500 CA WRITTEN UNKNOWN 00000600 C REVISED RDK; DEC 1980; IMPLEMENTED IN SPARC WITH NO CHANGES. 00000700 CA 00000800 CA 00000900 CA CALL MPOLRT (XCOF, COF, M, ROOTR, ROOTI, IER) 00001200 CA 00001300 CA 00001305 CA IN/OUT ARG TYPE DESCRIPTION 00001310 CA 00001320 CA IN XCOF R4 VECTOR OF M+1 COEFFICIENTS OF THE POLY- 00001400 CA NOMIAL ORDERED FROM SMALLEST TO LARGEST 00001500 CA POWER 00001600 CA IN COF R4 WORKING VECTOR OF LENGTH M+1 00001700 CA IN M I4 ORDER OF POLYNOMIAL 00001800 CA OUT ROOTR R4 RESULTANT VECTOR OF LENGTH M CONTAINING 00001900 CA REAL ROOTS OF THE POLYNOMIAL 00002000 CA OUT ROOTI R4 RESULTANT VECTOR OF LENGTH M CONTAINING 00002100 CA THE CORRESPONDING IMAGINARY ROOTS OF THE 00002200 CA POLYNOMIAL 00002300 CA OUT IER I4 ERROR CODE WHERE 00002400 CA = 0, NO ERROR 00002500 CA = 1, M LESS THAN ONE 00002600 CA = 2, M GREATER THAN 36 00002700 CA = 3, UNABLE TO DETERMINE ROOT WITH 500 00002800 CA INTERATIONS ON 5 STARTING VALUES 00002900 CA = 4, HIGH ORDER COEFFICIENT IS ZERO 00003000 CA 00003100 CA 00003200 CA THIS SUBROUTINE COMPUTES THE REAL AND COMPLEX ROOTS OF A REAL 00003300 CA POLYNOMIAL. 00003400 CA 00003500 CA 00003600 CA REMARKS: 00003700 CA 00003800 CA LIMITED TO 36TH ORDER POLYNOMIAL OR LESS. FLOATING POINT OVER- 00003900 CA FLOW MAY OCCUR FOR HIGH ORDER POLYNOMIALS BUT WILL NOT AFFECT 00004000 CA THE ACCURACY OF THE RESULTS. 00004100 CA 00004200 CA 00004300 CA SUBROUTINES REQUIRED: NONE 00004400 CA 00004500 CA 00004600 CA METHOD: 00004700 CA 00004800 CA NEWTON-RAPHSON ITERATIVE TECHNIQUE. THE FINAL ITERATIONS 00004900 CA ON EACH ROOT ARE PERFORMED USING THE ORIGINAL POLYNOMIAL 00005000 CA RATHER THAN THE REDUCED POLYNOMIAL TO AVOID ACCUMULATED 00005100 CA ERRORS IN THE REDUCED POLYNOMIAL. 00005200 C 00005300 SUBROUTINE MPOLRT(XCOF,COF,M,ROOTR,ROOTI,IER) 00005400 DIMENSION XCOF(1),COF(1),ROOTR(1),ROOTI(1) 00005500 C 00005600 DOUBLE PRECISION XO,YO,X,Y,XPR,YPR,UX,UY,V,YT,XT,U,XT2,YT2, 00005700 * SUMSQ,DX,DY,TEMP,ALPHA 00005800 C 00005900 C ...............................................................00006000 C 00006100 C IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE 00006200 C C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION 00006300 C STATEMENT WHICH FOLLOWS. 00006400 C 00006500 C DOUBLE PRECISION XCOF,COF,ROOTR,ROOTI 00006600 C 00006700 C THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS 00006800 C APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS 00006900 C ROUTINE. 00007000 C THE DOUBLE PRECISION VERSION MAY BE MODIFIED BY CHANGING THE 00007100 C CONSTANT IN STATEMENT 78 TO 1.0D-12 AND IN STATEMENT 122 TO 00007200 C 1.0D-10. THIS WILL PROVIDE HIGHER PRECISION RESULTS AT THE 00007300 C COST OF EXECUTION TIME 00007400 C 00007500 C ...............................................................00007600 C 00007700 IFIT=0 00007800 N=M 00007900 IER=0 00008000 IF(XCOF(N+1)) 10 , 40 , 10 00008100 10 IF(N) 20 , 20 , 60 00008200 C 00008300 C SET ERROR CODE TO 1 00008400 C 00008500 20 IER=1 00008600 30 RETURN 00008700 C 00008800 C SET ERROR CODE TO 4 00008900 C 00009000 40 IER=4 00009100 GO TO 30 00009200 C 00009300 C SET ERROR CODE TO 2 00009400 C 00009500 50 IER=2 00009600 GO TO 30 00009700 60 IF(N-36) 70 , 70 , 50 00009800 70 NX=N 00009900 NXX=N+1 00010000 N2=1 00010100 KJ1 = N+1 00010200 C 00010300 DO 80 L=1,KJ1 00010400 MT=KJ1-L+1 00010500 80 COF(MT)=XCOF(L) 00010600 C 00010700 C SET INITIAL VALUES 00010800 C 00010900 90 XO=.00500101 00011000 YO=0.01000101 00011100 C 00011200 C ZERO INITIAL VALUE COUNTER 00011300 C 00011400 IN=0 00011500 100 X=XO 00011600 C 00011700 C INCREMENT INITIAL VALUES AND COUNTER 00011800 C 00011900 XO=-10.0*YO 00012000 YO=-10.0*X 00012100 C 00012200 C SET X AND Y TO CURRENT VALUE 00012300 C 00012400 X=XO 00012500 Y=YO 00012600 IN=IN+1 00012700 GO TO 120 00012800 110 IFIT=1 00012900 XPR=X 00013000 YPR=Y 00013100 C 00013200 C EVALUATE POLYNOMIAL AND DERIVATIVES 00013300 C 00013400 120 ICT=0 00013500 130 UX=0.0 00013600 UY=0.0 00013700 V =0.0 00013800 YT=0.0 00013900 XT=1.0 00014000 U=COF(N+1) 00014100 IF(U) 140 , 290 , 140 00014200 C 00014300 140 DO 150 I=1,N 00014400 L =N-I+1 00014500 TEMP=COF(L) 00014600 XT2=X*XT-Y*YT 00014700 YT2=X*YT+Y*XT 00014800 U=U+TEMP*XT2 00014900 V=V+TEMP*YT2 00015000 FI=I 00015100 UX=UX+FI*XT*TEMP 00015200 UY=UY-FI*YT*TEMP 00015300 XT=XT2 00015400 150 YT=YT2 00015500 C 00015600 SUMSQ=UX*UX+UY*UY 00015700 IF(SUMSQ) 160 , 240 , 160 00015800 160 DX=(V*UY-U*UX)/SUMSQ 00015900 X=X+DX 00016000 DY=-(U*UY+V*UX)/SUMSQ 00016100 Y=Y+DY 00016200 170 IF(DABS(DY)+DABS(DX)-1.0D-05) 220 , 180 , 180 00016300 C 00016400 C STEP ITERATION COUNTER 00016500 C 00016600 180 ICT=ICT+1 00016700 IF(ICT-500) 130 , 190 , 190 00016800 190 IF(IFIT) 220 , 200 , 220 00016900 200 IF(IN-5) 100 , 210 , 210 00017000 C 00017100 C SET ERROR CODE TO 3 00017200 C 00017300 210 IER=3 00017400 GO TO 30 00017500 C 00017600 220 DO 230 L=1,NXX 00017700 MT=KJ1-L+1 00017800 TEMP=XCOF(MT) 00017900 XCOF(MT)=COF(L) 00018000 230 COF(L)=TEMP 00018100 C 00018200 ITEMP=N 00018300 N=NX 00018400 NX=ITEMP 00018500 IF(IFIT) 260 , 110 , 260 00018600 240 IF(IFIT) 250 , 100 , 250 00018700 250 X=XPR 00018800 Y=YPR 00018900 260 IFIT=0 00019000 270 IF(DABS(Y)-1.0D-4*DABS(X)) 300 , 280 , 280 00019100 280 ALPHA=X+X 00019200 SUMSQ=X*X+Y*Y 00019300 N=N-2 00019400 GO TO 310 00019500 290 X=0.0 00019600 NX=NX-1 00019700 NXX=NXX-1 00019800 300 Y=0.0 00019900 SUMSQ=0.0 00020000 ALPHA=X 00020100 N=N-1 00020200 310 COF(2)=COF(2)+ALPHA*COF(1) 00020300 C 00020400 320 DO 330 L=2,N 00020500 330 COF(L+1)=COF(L+1)+ALPHA*COF(L)-SUMSQ*COF(L-1) 00020600 C 00020700 340 ROOTI(N2)=Y 00020800 ROOTR(N2)=X 00020900 N2=N2+1 00021000 IF(SUMSQ) 350 , 360 , 350 00021100 350 Y=-Y 00021200 SUMSQ=0.0 00021300 GO TO 340 00021400 360 IF(N) 30 , 30 , 90 00021500 END 00021600