mblakele@jarthur.claremont.edu (Tad Blakeley) (02/22/91)
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
ENDwscott@en.ecn.purdue.edu (Wayne H Scott) (02/22/91)
I havn't posted this in a while, and some people seem to have
missed it so...
Article 1164 of comp.sys.handhelds:
Newsgroups: comp.sys.handhelds
Path: en.ecn.purdue.edu!wscott
From: wscott@ecn.purdue.edu (Wayne H Scott)
Subject: polynomial routines ver 3.1
Message-ID: <1990Sep3.153439.4558@ecn.purdue.edu>
Organization: Purdue University Engineering Computer Network
Date: Mon, 3 Sep 90 15:34:39 GMT
Here it is, my polynomial routines version 3.2
Changes from version 3.1:
- Faster PMUL
- RT now works with complex cooeffients
- various bug fixes
These routines are in the public domain, but I ask that if you use any of
them in one of your programs you give me credit. I am also
not responsible for any damage caused by these programs.
This package include the following programs.
TRIM Strip leading zeros from polynomial object.
IRT Invert root program. Given n roots it return a nth degree polynomial.
PDER Derivative of a polynomial.
RDER Derivative of a rational function.
PF Partial Fractions. (Handles multiple roots!)
FCTP Factor polynomial
RT Find roots of any polynomial
L2A Convert a list to an array and back.
PADD Add two polynomials
PMUL Multiply two polynomials.
PDIV Divide two polynomials.
EVPLY Evalulate a polynomial at a point.
COEF Given an equation return a polynomial list.
These programs are for the HP-48sx. I have a version of them that works
correctly on the HP-28. Send me mail if you want it.
I think people will find these very useful and work as I say, but if you
find any bugs please send me E-Mail. Comments are also welcome.
Some of these routines could be faster (PF, PMUL, ...) tell me if you know
how to speed them up.
_______________________________________________________________________________
Wayne Scott | INTERNET: wscott@en.ecn.purdue.edu
Electrical Engineering | BITNET: wscott%ea.ecn.purdue.edu@purccvm
Purdue University | UUCP: {purdue, pur-ee}!en.ecn.purdue.edu!wscott
_______________________________________________________________________________
These programs all work on polynomials in the follows form:
3*X^3-7*X^2+5 is entered as { 3 -7 0 5 }
Reasons why I use lists instead of arrays include:
* lists look better on the stack. (multi-line)
* The interactive stack can quickly put serveral items in a
list.
* Hitting EVAL will quickly return the elements on a list.
* the '+' key will add a element to the end or beginning of a list
or concat two lists
* Internally the main program that needs to be fast (BAIRS)
does 100% stack maniplations so speed of arrays vs. lists was
not a major issue.
* It would be easier for later releases to include symbolic
polynomials.
so going down the list...
The first program is FCTP. (factor polynomial)
When it is passed the cooeficients of a polynomial in a list it returns the
factor of that polynomal. ex:
1: { 1 -17.8 99.41 -261.218 352.611 -134.106 }
FCTP
3: { 1 -4.2 2.1 }
2: { 1 -3.3 6.2 }
1: { 1 -10.3 }
This tells us that X^5-17.8*X^4+99.41*X^3-261.218*X^2+352.611*X-134.106
factors to (X^2-4.2*X+2.1)*(X^2-3.3*X+6.2)*(X-10.3)
Neat!
The next program is RT. (Roots)
If given a polynmoial it return its roots. ex:
1: { 1 -17.8 99.41 -261.218 352.611 -134.106 }
RT
5: 3.61986841536
4: .58013158464
3: (1.65, 1.8648056199)
2: (1.65, -1.8648056199)
1: 10.3
Very Useful!
RT with work with complex cooeffients in the polynomial.
These programs use the BAIRS program which performs Bairstow's method of
quadratic factors and QUD with does the quadratic equation.
TRIM is used to strip the leading zeros from a polynomial list.
{0 0 3 0 7 } TRIM => { 3 0 7 }
RDER will give the derivative of a rational function.
ie. d x + 4 -X^2 - 8*x + 31
-- ------------- = --------------------------------
dx x^2 - 7*x + 3 x^4 - 14*x^3 + 55*x^2 - 42*x + 9
2: { 1 4 }
1: { 1 -7 3 }
RDER
2: { -1 -8 31 }
1: { 1 -14 55 -42 9 }
I don't know about you but I think it's useful.
IRT will return a polynomial whose roots you specify.
ie. If a transfer function has zeros at -1, 1 and 5 the function
is x^3 - 5*x^2 - x + 5
1: { -1 1 5 }
IRT
1: { 1 -5 -1 5 }
PDER will return the derivtive of a polynomial.
.ie The d/dx (x^3 - 5*x^2 - x + 5) = 3*x^2 - 10*x - 1
1: { 1 -5 -1 5 }
PDER
1: { 3 -10 -1 }
PF will do a partial fraction expansion on a transfer function.
.ie s + 5 1/18 5/270 2/3 1/9 2/27
----------------- = ----- + ----- - ------- - ------- - -----
(s-4)(s+2)(s-1)^3 (s-4) (s+2) (s-1)^3 (s-1)^2 (s-1)
2: { 1 5 }
1: { 4 -2 1 1 1 }
PF
1: { 5.5555e-2 1.85185e-2 -.6666 -.11111 -.074074 }
This program expects the polynomial of the numerator to be on level 2 and
a list with the poles to be on level 1. Repeated poles are suported but
they must be listed in order. The output is a list of the values of the
fraction in the same order as the poles were entered.
PADD, PMUL, and PDIV are all obvious, they take two polynomial lists off
the stack and perform the operation on them.
PDIV returns the polynomial and the remainder polynomial.
L2A converts a list to and array. (and back)
1: { 1 2 3 }
L2A
1: [ 1 2 3 ]
L2A
1: { 1 2 3 }
EVPLY evalutates and polynomial at a point.
x^3 - 3*x^2 +10*x - 5 | x=2.5 = 16.875
2: { 1 -3 10 -5 }
1: 2.5
EVPLY
1: 16.875
P.S. Many thanks to Mattias Dahl & Henrik Johansson for fixs they have
made.
%%HP: T(3)A(R)F(.);
DIR
PDIV
\<< DUP SIZE 3 ROLLD OBJ\-> \->ARRY SWAP OBJ\-> \->ARRY \-> c b a
\<< a b
IF c 1 SAME THEN
OBJ\-> DROP / OBJ\-> 1 GET \->LIST { 0 }
ELSE
WHILE
OVER SIZE 1 GET c \>=
REPEAT DIVV
END
DROP \-> d
\<< a SIZE 1 GET c 1 - - IF DUP NOT THEN
1 END \->LIST d OBJ\-> OBJ\-> DROP \->LIST
\>>
END
\>>
\>>
TRIM
\<< OBJ\-> \-> n
\<< n
WHILE ROLL DUP ABS NOT n 1 - AND
REPEAT DROP 'n' DECR
END n ROLLD
n \->LIST
\>>
\>>
RDER
\<< \-> F G
\<< G F
PDER PMUL G PDER {
-1 } PMUL F PMUL
PADD G G PMUL
\>>
\>>
IRT
\<< OBJ\-> \-> n
\<<
IF n 0
>
THEN 1
n
START
n ROLL { 1 } SWAP
NEG +
NEXT
ELSE {
1 }
END
IF n 1
>
THEN 2
n
START
PMUL
NEXT
END
\>>
\>>
PDER
\<< OBJ\-> \-> n
\<< 1 n
FOR i n
ROLL n i - *
NEXT
DROP
IF n 1
==
THEN {
0 }
ELSE n
1 - \->LIST
END
\>>
\>>
PF
\<< MAXR { }
\-> Z P OLD LAST
\<< 1 P
SIZE
FOR I P
I GET \-> p1
\<<
IF p1 OLD \=/
THEN Z p1 EVPLY 1 P
SIZE
FOR J
IF P J GET P I
GET \=/
THEN p1 P J GET
- /
END
NEXT p1 'OLD' STO
{ } 'LAST' STO
ELSE
IF { } LAST SAME
THEN 1 { } 1 P
SIZE
FOR K P K GET
IF DUP p1 ==
THEN DROP
ELSE +
END
NEXT IRT Z SWAP
ELSE LAST OBJ\->
DROP
END RDER DUP2 5
PICK 1 + 3 ROLLD 3
\->LIST 'LAST' STO p1
EVPLY SWAP p1 EVPLY
SWAP / SWAP ! /
END
\>>
NEXT P
SIZE \->LIST
\>>
\>>
FCTP
\<<
IF DUP
SIZE 3 >
THEN DUP
BAIRS SWAP OVER
PDIV DROP FCTP
END
\>>
RT
\<< TRIM DUP
SIZE \-> n
\<<
IF n 3
>
THEN
DUP BAIRS SWAP OVER
PDIV DROP \-> A B
\<< A
RT B RT
\>>
ELSE
IF n
2 >
THEN
QUD
ELSE
LIST\-> DROP NEG SWAP
/
END
END
\>>
\>>
L\178A
\<<
IF DUP
TYPE 5 ==
THEN OBJ\->
\->ARRY
ELSE OBJ\->
1 GET \->LIST
END
\>>
PADD
\<< DUP2 SIZE
SWAP SIZE \-> A B nB
nA
\<< A L\178A B
L\178A
IF nA
nB <
THEN
SWAP
END
IF nA
nB \=/
THEN 1
nA nB - ABS
START
0
NEXT
END nA
nB - ABS 1 + ROLL
OBJ\-> 1 GET nA nB -
ABS + \->ARRY + L\178A
\>>
\>>
PMUL
\<< DUP2 SIZE
SWAP SIZE \-> X Y ny
nx
\<<
1 nx ny + 1 - FOR I
0
NEXT
1 nx FOR I
1 ny FOR J
I J + 1 - ROLL
X I GET Y J GET * +
I J + 1 - ROLLD
NEXT
NEXT
{ }
1 nx ny + 1 - START
SWAP +
NEXT
\>>
\>>
EVPLY
\<< OVER
IF DUP
TYPE 5 ==
THEN SIZE
ELSE SIZE
1 GET
END \-> a r
n
\<< a 1 GET
IF n 1
>
THEN 2
n
FOR i
r * a i GET +
NEXT
END
\>>
\>>
COEF
\<< \-> E n
\<< 0 n
FOR I 0
'X' STO E EVAL 'X'
PURGE E 'X' \.d 'E'
STO I ! /
NEXT 2
n 1 +
FOR I I
ROLL
NEXT n
1 + \->LIST
\>>
\>>
EQ 1
DIVV
\<< DUP 1 GET \-> a
b c
\<< a 1 GET c /
DUP b * a SIZE RDM
a SWAP - OBJ\-> 1
GETI 1 - PUT \->ARRY
SWAP DROP b
\>>
\>>
QUD
\<< LIST\-> \->ARRY
DUP 1 GET / ARRY\->
DROP ROT DROP SWAP
2 / NEG DUP SQ ROT
- \v/ DUP2 + 3 ROLLD
-
\>>
BAIRS
\<< OBJ\-> 1 1 \-> n
R S
\<<
DO 0 n 1 +
PICK 0 0 0 4 PICK 5
n + 7
FOR J J
PICK R 7 PICK * + S
8 PICK * + 7 ROLL
DROP DUP 6 ROLLD R
3 PICK * + S 4 PICK
* + 5 ROLL DROP -1
STEP 3
PICK SQ 3 PICK 6
PICK * -
IF DUP 0
==
THEN DROP
1 1
ELSE 6
PICK 6 PICK * 5
PICK 9 PICK * -
OVER / 4 PICK 9
PICK * 8 PICK 7
PICK * - ROT /
END DUP
'S' STO+ SWAP DUP
'R' STO+
UNTIL (0,1) * +
ABS .000000001 < 7
ROLLD 6 DROPN
END n DROPN
1 R NEG S NEG 3
\->LIST
\>>
\>>
END
--
_______________________________________________________________________________
Wayne Scott | INTERNET: wscott@en.ecn.purdue.edu
Electrical Engineering | BITNET: wscott%ea.ecn.purdue.edu@purccvm
Purdue University | UUCP: {purdue, pur-ee}!en.ecn.purdue.edu!wscott
--
_________________________________________________________________________
Wayne Scott | INTERNET: wscott@ecn.purdue.edu
Electrical Engineering | BITNET: wscott%ecn.purdue.edu@purccvm
Purdue University | UUCP: purdue, pur-ee}!ecn.purdue.edu!wscott