Generating a Polynomial Given Its Roots

[Eddie W. Shore]

Generate the coefficients of a polynomial (up to the order 4) with the roots a_0, a_1, a_2, and a_3. The resulting polynomial is:

p(x) = (x - a_0) * (x - a_1) * (x - a_2) * (x - a_3) * (x - a_4)

p(x) = r_4 * x^4 + r_5 * x^3 + r_6 * x^2 + r_7 * x + r_8

The default is a polynomial where the lead coefficient is positive. If you want a polynomial where the lead coefficient is negative, multiply every coefficient by -1.

Instructions

Store the four roots in registers R00, R01, R02, and R03 respectively. Run POLY4. Coefficients are shown briefly as they are calculated. They are can be recalled by the registers in decreasing order of x: R04, R05, R06, R07, and R08.

DM 41L and HP 41C Program: POLY4

Code:


01  LBL^T POLY4
02 1
03 STO 04
04 PSE 
05 RCL 00
06 CHS
07 RCL 01
08 -
09 RCL 02
10 -
11 RCL 03
12 -
13 STO 05
14 PSE
15 RCL 01
16 RCL 02
17 +
18 RCL 03
19 +
20 RCL 00
21 * 
22 RCL 02
23 RCL 03
24 +
25 RCL 01
26 *
27 +
28 RCL 02
29 RCL 03
30 *
31 +
32 STO 06
33 PSE
34 RCL 01
35 RCL 02
36 *
37 RCL 01
38 RCL 03
39 *
40 +
41 RCL 02
42 RCL 03
43 *
44 +
45 RCL 00
46 *
47 CHS
48 RCL 01
49 RCL 02
50 *
51 RCL 03
52 *
53 -
54 STO 07
55 PSE
56 RCL 00
57 RCL 01
58 *
59 RCL 02
60 *
61 RCL 03
62 *
63 STO 08
64 RTN

Example

Roots x = -3, x = 3, x= 4, and x= 6

Coefficients:
R04 = 1
R05 = -10
R06 = 15
R07 = 90
R08 = -216

Polynomial: p(x) = x^4 - 10 * x^3 + 15 * x^2 + 90 * x - 216

Blog post: https://edspi31415.blogspot.com/2018/12/...omial.html

[Namir]

Your post and code is very interesting!! I am working on a version that will build polynomials up to order of 19. I have tested the code using Excel VBA. Next is coding the program for HP-41CX. I may also write a version for HP-71B.

Namir

[Thomas Klemm]

If you really only want to do that up to order 4 you can use:

Code:

LBL "POLY4"
CHS   ENTER↑   X<> 00   ST+ 00
RCL Y   *   X<> 01   ST+ 01
RCL Y   *   X<> 02   ST+ 02
RCL Y   *   X<> 03   ST+ 03
END

Example:

CLRG
-3 XEQ "POLY4"
3 R/S
4 R/S
6 R/S

The coefficients of the polynomial can then be found in the registers:

R 00:  -10.0000
R 01:   15.0000
R 02:   90.0000
R 03: -216.0000

---

But we can do better than loop unrolling and use for the general case:

Code:

01 LBL "POLY"
02 CHS
03 ENTER↑
04 GTO 01
05 LBL 00
06 RCL Y
07 *
08 LBL 01
09 X<> IND Z
10 ST+ IND Z
11 ISG Z
12 GTO 00
13 RDN
14 RDN
15 FRC
16 1 E-3
17 +
18 END

Example:

CLRG
CLST
-3 XEQ "POLY"
3 R/S
4 R/S
6 R/S

The coefficients of the polynomial can again be found in the registers:

R 00:  -10.0000
R 01:   15.0000
R 02:   90.0000
R 03: -216.0000


In both cases the leading coefficient is always 1.

Cheers
Thomas

[Namir]

Here is my version of the HP-41C program that can handle up to 19 roots.

The code basically implements repeatedly multiplying a polynomial with a simple 1st degree polynomial. The initial form of the polynomial is also a simple first degree polynomial. As you enter more roots, the degree of the polynomial increases.

Code:


Pseudo-Code
========
A(1) = 1
A(2) = -first root
Do
  B = input("Enter root?")
  IF B = 0 Then Exit Loop
  B = -B
  For I=2 to N+1
    C(I) = A(I) + B * A(I-1)
  Next I
  C(N+2) = A(N+1)*B
  N = N + 1
  For I = 1 To N + 1
    A(I) = C(I)
  Next I
Loop

Display array A(1) to A(N+1)

Memory Map
==========

R00 = I loop control variable
R01 = A(1)
R02 = A(2)
...
R20 = A(20)
R21 = C(1)
R22 = C(2)
...
R40 = C(20)
R41 = N
R42 = B
R43 = index I

HP-41C Listing
==============

01 LBL "POLY"        
02 LBL A        
03 1        
04 STO 41        
05 STO 01        
06 STO 21        
07 FIRST ROOT?        
08 PROMPT        
09 CHS        
10 STO 02        
11 STO 22        
12 LBL 00    # Main loop
13 ROOT? 0=END        
14 PROMPT        
15 X=0?        
16 GTO 01    # go to display results
17 CHS        
18 STO 42    # store minus root
19 XEQ 09        
20 ISG 00    # store value for loop conttrol variable
21 STO X    # a virtual NOP
22 LBL 02    # stert inner loop
23 RCL 00        
24 INT        
25 STO 43    # store index I
26 1        
27 -        
28 RCL IND X    # Get A(I-1)    
29 RCL 42        # Get B
30 *        
31 RCL IND 43    # Get A(I)    
32 +        
33 RCL 00        
34 INT        
35 20        
36 +        
37 X<>Y        
38 STO IND Y    # Store in C(I)    
39 ISG 00        
40 GTO 02    # End of loop
41 RCL 41        
42 1        
43 +        
44 RCL IND X        
45 RCL 42        
46 *        
47 RCL 41        
48 2        
49 +        
50 X<>Y        
51 STO IND Y        
52 XEQ 09        
53 LBL 03    # Loop for A(I) = C(I)
54 RCL 00        
55 INT        
56 STO 43        
57 20        
58 +        
59 RCL IND X        
60 STO IND 43        
61 ISG 00        
62 GTO 03        
63 1        
64 STO+ 41    # N = N + 1 increment polynomial order
65 GTO 00        
66 LBL 01        
67 XEQ 09        
68 LBL 04        
69 RCL 00        
70 INT        
71 RCL IND X        
72 R/S        
73 ISG 00        
74 GTO 04        
75 END        
76 PROMPT        
77 RTN        
78 LBL 09    # set up loop control variable
79 RCL 41        
80 1        
81 +        
82 1.00E+03        
83 /        
84 1        
85 +        
86 STO 00        
87 RTN

Run the program by executing XEQ "POLY" (you can instead press the key [A] when you are inside the program POLY):

1) The first prompt will ask you for the first root. Enter a root.
2) Subsequent prompts will ask you for additional roots OR to enter 0 when you are done.
3) The program displays the coefficients of the polynomial starting with the coefficient of the highest term. Press [R/S] to view the next coefficient of the sequentially lower term. When you have viewed all of the coefficients, the program displays END.

Inspect the value in register 41 to view the degree of the resulting polynomial.

Enjoy!

[Thomas Klemm]

Here's the program adapted for the HP-11C:

Code:

001-42,21,11    LBL A
002-      16    CHS
003-      36    ENTER
004-    22 1    GTO 1
005-42,21, 0    LBL 0
006-      34    x<>y
007-      20    ×
008-   43 36    LSTx
009-      34    x<>y
010-42,21, 1    LBL 1
011-   42 23    x<> (i)
012-44,40,24    STO+ (i)
013-    42 6    ISG
014-    22 0    GTO 0
015-   45 25    RCL I
016-   42 44    FRAC
017-      26    EEX
018-       3    3
019-      16    CHS
020-      40    +
021-   44 25    STO I
022-   43 32    RTN

Example:

CLEAR REG
CLx
STO I
-3 A
3 R/S
4 R/S
6 R/S

The coefficients of the polynomial can be found in the registers:

R 0:  -10.0000
R 1:   15.0000
R 2:   90.0000
R 3: -216.0000

Here as well the leading coefficient is always 1.

Cheers
Thomas

[Eddie W. Shore]

(12-22-2018 11:44 PM)Namir Wrote:  Your post and code is very interesting!! I am working on a version that will build polynomials up to order of 19. I have tested the code using Excel VBA. Next is coding the program for HP-41CX. I may also write a version for HP-71B.

Namir

Thank you! I love your work.

[Namir]

You have a fantastic web site too!!!

Namir

[Namir]

I added pseudo-code to my original post.

Namir

[Namir]

As I promised, here is the HP-71B version of the program. The DATA statement contain the number of roots, followed by the list of roots. This data can appear on one or more DATA statements. Run the program to view the coefficients along with their associated power of X:

Code:

10 REM POLY COEFFS FROM ROOTS
20 DIM A(101),C(101)
25 REM DATA NUMBER_OF_ROOTS
30 DATA 3
35 REM DATA STMT HAS LIST OF ROOTS
40 DATA 1,2,3
50 READ M, B
60 A(1) = 1@ C(1) = 1
70 A(2) = -B@ C(2) = -B
80 N = 1
90 FOR K=2 TO M
100 READ B @ B=-B
110 FOR I = 2 TO N + 1
120 C(I) = A(I) + B * A(I - 1)
130 NEXT I
140 C(N + 2) = A(N + 1) * B
150 N = N + 1
160 FOR I = 1 TO N + 1
170 A(I) = C(I)
180 NEXT I
190 NEXT K
195 M = N
200 FOR I=1 TO N
210 DISP A(I);"* X^";M @ WAIT 3
215 M = M - 1
220 NEXT I
230 DISP A(N+1)
240 END

[Thomas Klemm]

(12-26-2018 12:09 AM)Namir Wrote:  As I promised, here is the HP-71B version of the program.

I was a bit surprised that the HP-71B is missing the PCOEF command of the HP-48G:

[-3 3 4 6]
PCOEF

[ 1 -10 15 90 -216 ]

Quote: The DATA statement contain the number of roots, followed by the list of roots.

So we have to modify the program for every new set of data?
That's not very practical.

Cheers
Thomas

[Namir]

(12-27-2018 10:35 PM)Thomas Klemm Wrote:  

(12-26-2018 12:09 AM)Namir Wrote:  As I promised, here is the HP-71B version of the program.

I was a bit surprised that the HP-71B is missing the PCOEF command of the HP-48G:

[-3 3 4 6]
PCOEF

[ 1 -10 15 90 -216 ]

Quote: The DATA statement contain the number of roots, followed by the list of roots.

The DATA statements provide the input. It's easier to change one or two numbers, add a number (then RUN the program) than to enter the whol list of roots again.

Namir

So we have to modify the program for every new set of data?
That's not very practical.

Cheers
Thomas