Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!usc!jarthur!jarthur.claremont.edu
From: mblakele@jarthur.claremont.edu (Tad Blakeley)
Newsgroups: comp.sys.handhelds
Subject: Root-finder for HP48
Message-ID: <10917@jarthur.Claremont.EDU>
Date: 21 Feb 91 21:01:45 GMT
Sender: mblakele@jarthur.Claremont.EDU
Organization: Harvey Mudd College, Claremont, CA 91711
Lines: 244


	Here's a beta version of PROOT, a root-finder for real-coefficient
polynomials.  Inputs are N, the degree of the polynomial, and A, its
coefficient matrix.  Example:  to find the roots of x^3+2x^2+3x+4,
enter [1 2 3 4] 'A' STO 3 'N' STO PROOT
	The routines were ported from a Fortran program, so I'm still
optimizing.  Eventually, if I can get the speed up to speed, I'll build
a root-locus plotter around PROOT.
	If you find this program useful, please try to optimize it.  The
number above the returned roots is a time for the PROOT call; it takes
about 16 seconds on my HP for a cubic polynomial.  If you can improve the
speed, accuracy, or size of PROOT, please let me know.  I'll post
improved versions as I get them written.
					-- tad
------------------------clip 'n' save-------------------------------------
%%HP: T(1)A(D)F(.);
DIR
  A
[ 1 2 3 4 ]
  N 3
  PROOT
    + TIME INIT
PART3 PART4
    ;
  TRIM
    + TIME - U V
(0,1) * + OBJ
 DROP
DROP KILL
    ;
  U
[ -.174685404255 -.174685404255 -1.65062919149 0 ]
  V
[ -1.54686888726 1.54686888726 0 0 ]
  INIT
    + N 1 + 1 
LIST
0 CON DUP 'V' STO
'U' STO N 1 + 'NC'
STO NC 1 + 1 
LIST
0 CON DUP 'C' STO
'B' STO -1 'IRREV'
STO A 'H' STO 0 0 0
'P' STO 'Q' STO 'R'
STO
    ;
  PP .174685404255
  QP 2.42331834482
  F -2.39280335436
  NL 3
  M 2
  NP 3
  E .0000000001
  PART70
    + 'NC' 2 STO-
      IF IRREV 0 <
      THEN Q INV
'QP' STO P Q 2 * /
'PP' STO
      ELSE Q 'QP'
STO P 2 / 'PP' STO
      END PP SQ QP
- 'F' STO
      IF F 0 <
      THEN U NC 1 +
PP NEG PUT NC PP
NEG PUT 'U' STO V
NC 1 + F NEG  PUT
DUP NC 1 + GET NEG
NC SWAP PUT 'V' STO
      ELSE U NC 1 +
PP PP ABS / PP ABS
F  + * NEG PUT DUP
NC 1 + GET QP SWAP
/ NC SWAP PUT 'U'
STO V NC 1 + 0 PUT
NC 0 PUT 'V' STO
      END 1 NC
      FOR I H I B I
2 + GET PUT 'H' STO
      NEXT PART4
    ;
  PART50
    + 'NC' 1 STO- V
NC 0 PUT 'V' STO
      IF IRREV 0 <
      THEN U NC R
INV PUT 'U' STO
      ELSE U NC R
PUT 'U' STO
      END 1 NC
      FOR I H I B I
1 + GET PUT 'H' STO
      NEXT PART4
    ;
  PART20
    + 1 1000
      FOR J NC 1 -
NC NP -
        FOR I B I H
I GET R B I 1 + GET
* + PUT 'B' STO C I
B I GET R C I 1 +
GET * + PUT 'C' STO
-1
        STEP
        IF B 1 GET
H 1 GET / ABS E 	
        THEN PART50
        END
        IF C 2 GET
0 ==
        THEN R 1 +
'R' STO
        ELSE R B 1
GET C 2 GET / - 'R'
STO
        END NC 1 -
NC NP -
        FOR I B I H
I GET P B I 1 + GET
* - Q B I 2 + GET *
- PUT 'B' STO C I B
I GET P C I 1 + GET
* - Q C I 2 + GET *
- PUT 'C' STO -1
        STEP 1 'T'
STO
        IF H 2 GET
0 
        THEN 1 'T'
STO+
        END
        IF B 2 GET
H T GET / ABS E 	 B
1 GET H 1 GET / ABS
E 	 AND
        THEN PART70
        END C 2 GET
B 2 GET - 'CBAR'
STO C 3 GET SQ CBAR
C 4 GET * - 'D' STO
        IF D 0 ==
        THEN P 2 -
'P' STO Q Q 1 + *
'Q' STO
        ELSE P B 2
GET C 3 GET * B 1
GET C 4 GET * - D /
+ 'P' STO
        END
      NEXT Q B 2
GET CBAR * B 1 GET
C 3 GET * - D / -
'Q' STO E 10 * 'E'
STO PART20
    ;
  PART19
    + .0000000001
'E' STO B NC H NC
GET PUT 'B' STO C
NC H NC GET PUT 'C'
STO B NC 1 + 0 PUT
'B' STO C NC 1 + 0
PUT 'C' STO NC 1 -
'NP' STO PART20
    ;
  PART4
    +
      IF NC 1 - 0
==
      THEN TRIM
      END
      IF NC 2 - 0
==
      THEN H 1 GET
H 2 GET / NEG 'R'
STO PART50
      END
      IF NC 3 - 0
==
      THEN H 2 GET
H 3 GET / 'P' STO H
1 GET H 3 GET / 'Q'
STO PART70
      END
      IF H NC 1 -
GET H NC GET / ABS
H 2 GET H 1 GET /
ABS - 0 

      THEN PART19
      END IRREV NEG
'IRREV' STO NC 2 /
'M' STO 1 M
      FOR I NC 1 +
I - 'NL' STO H NL
GET 'F' STO H NL H
I GET PUT 'H' STO H
I F PUT 'H' STO
      NEXT
      IF Q 0 
      THEN P Q /
'P' STO Q INV 'Q'
STO
      ELSE 0 'P'
STO
      END
      IF R 0 
      THEN R INV
'R' STO
      END PART19
    ;
  PART3
    + H 1 GET
      IF 0 ==
      THEN 'NC' 1
STO- V NC 0 PUT 'V'
STO U NC 0 PUT 'U'
STO 1 NC
        FOR I H I H
I 1 + GET PUT 'H'
STO
        NEXT PART3
      END
    ;
  R -1.65062919149
  Q 2.42331834482
  P .34937080851
  H
[ .34937080851 .34937080851 1 1 ]
  IRREV 1
  B
[ -.00000000023 2.42331834482 .34937080851 1 0 ]
  C
[ -7.54537830759 4.57121341744 -1.30125838298 1 0 ]
  NC 1
  D 3.43041707098
  CBAR
-.17701875726
  T 2
  I 1
END