Polynomial Ops Program

[Namir]

Polynomial Operations

This program supports the following polynomial operations:

1) Input and viewing of polynomial coefficients.
2) Calculating polynomial value at scalar X.
3) Calculating polynomial's first derivative at scalar X.
4) Calculating polynomial's second derivative at scalar X.
5) Calculating definite integrals for polynomials.
6) Calculating real polynomial roots using Newton's method.
7) Deflating a polynomial.
8) REplacing the original polynomial with its deflated version.

I included a ZIP file with a .raw file that contains the programs. You can load it into an HP-41C emulator that accepts .raw files.

Synthetic Division

P(x) = a0 + a1*X + a2*X^2 + ... + an*X^N
Q(x) = (x-X0)(b0 + b1*X + b2*X^2 + ... + b<n-1>*X^N-1) + R
Where:
b<n-1> = a<n>
b<k> = b<k+1> * X0 + a<k+1> for k = N-2, N-3,...,1,0
R = b0*X0 + a0

Polynomial Derivatives

P'(X) = a1 + 2*a2*X + 3*a3*X^2 + ... + N*a<n>*X^N-1
P''(X) = 2*a2 + 6*a3*X + 12*a4*X^2 + ... + N*(N-1)*a<n>*X^N-2

Polynomial Integrals

T(X) = Integeral(P(X),A,B) = a0*(B-A) + a1*(B^2-A^2)/2 + a2*(B^3-A^3)/3 + ... + a<n>*(B^(N+1)-A^(N+1))/(N+1)

Real Polynomial Roots

Root of polynomial uses Newton's method:
Xnew = X - P(X) / P'(X)


Memory Map

R00 = N
R01 = X0, counter
R02 = index for a() = 10 + (n+10)/1000 or n+10 + 10/1000
R03 = index for b() = n+11 + (2n+10)/1000 or 2n+10 + (n+11)/1000
R04 = X, A
R05 = Power of X, B
R06 = Sum
R07 = Power
R08 = counter, X
R09 = tolerance
R10 = a(0)
R11 = a(1)
...
Rn+10 = a(n)
Rn+11 = b(0)
Rn+12 = b(1)
Rn+13 = b(2)
...
Rn-1+n+11 = R2n+10 = b(n-1)

Flag 0 = Set means deflated polynomial exist.
Flag 01 - used by root-seeking routine to calculate P(X)
Flag 02 - used

Listing

Code:

1 LBL "POLY1"
2 LBL A        # input OR view polynomial coefficients
3 0
4 STO 01        # store 0 in counter
5 CF 29
6 FIX 0
7 "N="
8 ARCL 00
9 SF 29
10 FIX 5
11 CF 22
12 PROMPT        # prompt for polynomial order
13 FS? 22        # overwrite value with user's input
14 STO 00
15 FS? 22        # clear flag 00 since user altered polynomial?
16 CF 00
17 XEQ 00        # store index to access a()
18 LBL 05        #-------------------  Start loop to prompt and enter a(i), i=0,1,...,n
19 "A<"
20 CF 29
21 FIX 0
22 ARCL 01
23 SF 29
24 FIX 5
25 |-">="
26 ARCL IND 02
27 CF 22
28 PROMPT
29 FS? 22        # overwrite value with user's input
30 STO IND 02
31 FS? 22        # clear flag 00 since user altered polynomial?
32 CF 00
33 ISG 01        # increment counter
34 STO X
35 ISG 02
36 GTO 05        # ------------------- end of input loop
37 RTN
38 LBL a        # Calculate a root for P(X)
39 "X?"
40 PROMPT
41 STO 08        # store X
42 CF 22
43 "TOLER?"
44 PROMPT
45 FC?22
46 1E-9        # default tolerance value
47 STO 03        # store tolerance value
48 0
49 STO 01        # initialize iteration counter
50 LBL 11        # ------------------- start of loop
51 1
52 STO+ 01        # increment iter counter
53 RCL 01
54 50    
55 X<=Y?        # Is 50 <= iter counter?
56 GTO 12        # exit loop
57 RCL 08
58 PSE
59 SF 01
60 XEQ c        # Calculate P'(X)
61 X=0?        # is P'(X)=0?
62 GTO 13        # exit loop
63 STO 09
64 RCL 08
65 SF 01
66 XEQ C        # Calculate P(X)
67 RCL 09
68 /
69 STO- 08        # X = X - P(X)/P'(X)
70 ABS
71 RCL 03
72 X<=Y?
73 GTO 11        # ------------------- end of loop
74 "RT="    
75 RCL 08
76 ARCL X
77 PROMPT        # display root
78 RTN
79 LBL 12
80 "MAX ITERS"
81 AVIEW
82 PSE
83 "AT X="
84 ARCL 04
85 PTOMPT
86 RTN
87 LBL 13
88 "NO ROOT"
89 PROMPT
90 RTN

91 LBL B        # deflate P(X)
92 "X0?"
93 PROMPT
94 STO 01
95 XEQ 01
96 XEQ 03
97 RCL 00
98 1
99 -
100 STO 08        # initialize index for b()
101 RCL IND 02
102 STO IND 03    # b(n-1) = a(n)
103 "B<"
104 CF 29
105 FIX 0
106 ARCL 08
107 SF 29
108 FIX 5
109 ":->="
110 ARCL IND 03
111 PROMPT        # display b(n-1)
112 1
113 STO- 02
114 LBL 10        # ------------------- Start loop for deflation
115 1
116 STO- 08        # decrement index for b()
117 RCL IND 03
118 RCL 01
119 *
120 RCL IND 02
121 +        # calculate b(k)
122 DSE 03
123 STO X
124 STO IND 03    # store b(k)
125 "B<"
126 CF 29
127 FIX 0
128 ARCL 08
129 SF 29
130 FIX 5
131 ":->="
132 ARCL IND 03
133 PROMPT        # display b(k)
134 DSE 02
135 GTO 10        # ------------------- end deflation loop
136 SF 00
137 RCL IND 03
138 RCL 01
139 *
140 RCL IND 02
141 +        # calculate residual R
142 "R="
143 ARCL X
144 PROMPT
145 RTN
146 LBL C        # Evaluate P(X)
147 XEQ 01        # Setup the index to access a()
148 "X?"
149 FC? 01        # prompt for X if flag 01 is clear
150 PROMPT
151 FS? 01
152 RCL 08        # provide X from caller routine
153 ENTER
154 ENTER
155 ENTER
156 SF 02
157 LBL 06        # --------------- start loop
158 RCL IND 02
159 FC?C 02
160 +
161 *
162 DSE 02
163 GTO 06        # --------------- end of loop
164 RCL IND 02
165 +
166 "P="
167 ARCL X
168 FC?C 01
169 PROMPT        # Display polynomial value if flag 01 is clear
170 RTN
171 LBL c        # Evaluate P'(X)
172 XEQ 00        # Setup the index to access a(0) to a(n)
173 1
174 +
175 STO 02        # Store adjusted index to skip accessing a(0)
176 "X?"
177 FC? 01
178 PROMPT
179 FS? 01
180 RCL 08        # provide X from caller routine
181 STO 04        # Store X
182 1
183 STO 05        # power of X = 1
184 STO 07
185 0
186 STO 06        # sum = 0
187 LBL 09        # ----------------------------------- start loop
188 RCL 07
189 RCL IND 02
190 *
191 RCL 05
192 *
193 STO+ 06        # sum = sum + power*a(i)Xpower
194 1
195 STO+ 07        # power = power + 1
196 RCL 04
197 STO* 05        # Xpower = X * Xpower
198 ISG 02        # increment coefficient accessing index
199 GTO 09        # ----------------------------------- end loop
200 RCL 06
201 "D1="
202 ARCL X
203 FC?C 01
204 PROMPT
205 RTN
206 LBL D        # Evaluate Q(X)
207 XEQ 03        # Setup the index to access b()
208 "X?"
209 PROMPT
210 ENTER
211 ENTER
212 ENTER
213 SF 02
214 LBL 07        # -------------------  start loop
215 RCL IND 03
216 FC?C 02
217 +
218 *
219 DSE 03
220 GTO 07        # -------------------  end loop
221 RCL IND 03
222 +
223 "Q="
224 ARCL X
225 PROMPT
226 RTN
227 LBL d        # Evaluate P''(X)
228 XEQ 00        # Setup the index to access a(0) to a(n)
229 2
230 +
231 STO 02        # Store adjusted index to skip accessing a(0) and a(1)
232 "X?"
233 PROMPT
234 STO 04        # Store X
235 1
236 STO 05        # power of X = 1
237 STO 07
238 0
239 STO 06        # sum = 0
240 LBL 19        # ----------------------------------- start loop
241 RCL 07
242 X^2
243 LastX
244 +
245 RCL IND 02
246 *
247 RCL 05
248 *
249 STO+ 06        # sum = sum + power*a(i)Xpower
250 1
251 STO+ 07        # power = power + 1
252 RCL 04
253 STO* 05        # Xpower = X * Xpower
254 ISG 02        # increment coefficient accessing index
255 GTO 19        # ----------------------------------- end loop
256 RCL 06
257 "D2="
258 ARCL X
259 PROMPT
260 RTN
261 LBL E        # Replace P(X) by Q(X)
262 FC? 00        # If flag 00 is clear, exit since the coefficients of Q(X) were not yet calculated
263 RTN
264 10
265 STO 02        # Set simple index to access a(i)
266 XEQ 02
267 LBL 08        # ------------------- start loop
268 RCL IND 03
269 STO IND 02    # b(i)=a(i)
270 1
271 STO+ 02        
272 ISG 03
273 GTO 08        # ------------------ end loop
274 CF 00
275 1
276 STO- 00        # n=n-1
277 RTN
278 LBL e        # Integrate P(X) between A and B
279 XEQ 00        # Setup the index to access a(0) to a(n)
280 "A^B?"
281 PROMPT
282 STO 04        # Store A
283 X<>Y
284 STO 05        # Store B
285 1
286 STO 07
287 0
288 STO 06        # sum = 0
289 LBL 29        # ----------------------------------- start loop
290 RCL 05
291 RCL 07
292 Y^X
293 RCL 04
294 RCL 07
295 Y^X
296 -
297 RCL 07
298 /
299 RCL IND 02
300 *
301 STO+ 06        # sum = sum + power*a(i)Xpower
302 1
303 STO+ 07        # power = power + 1
304 ISG 02        # increment coefficient accessing index
305 GTO 29        # ----------------------------------- end loop
306 RCL 06
307 "INTTG="
308 ARCL X
309 PROMPT
310 RTN
311 LBL 00        # Setup the index to access a() in ascending order (index= 10 + (n+10)/1000)
312 RCL 00
313 10
314 STO Z
315 +
316 1E3
317 /
318 +
319 STO 02
320 RTN
321 LBL 01        # Setup the index to access a() in descending order (index = n+10 + 10/1000)
322 RCL 00
323 10.01
324 +
325 STO 02
326 RTN
327 LBL 02        # Setup the index to access b() in ascending order (index =  n+11 + (2n+10)/1000)
328 RCL 00
329 STO+ X
330 10
331 +
332 1E 3
333 /
334 RCL 00
335 +
336 11
337 +
338 STO 03        
339 RTN
340 LBL 03        # store index to access b()  in descending order (index = 2n+10 + (n+11)/1000)
341 RCL 00
342 STO+ X
343 10
344 +
345 RCL 00
346 11
347 +
348 1E 3
349 /
350 +
351 STO 03
352 RTN

Usage

1) To enter the coefficients a(), press [A]. The program displays the polynomial order and then the value of the coefficients. To change any of these values, enter a new value. Otherwise, press [R/S] to view/edit the next value.
2) To deflate the polynomial P(x), press {B}. The program will prompt you for X0 and then display the coefficients b() of the deflated polynomial Q(x).
3) To calculate P(X), press [C]. The program prompt you for X and then calculate P(X).
4) To calculate P'(X) press [f][C]. The program prompt you for X and then calculate P'(X).
5) To calculate P''(X) press [f][D]. The program prompt you for X and then calculate P''(X).
6) To calculate the integral of P(X) press [f][E]. The program prompt you for integral range [A, B] and then calculate the integral of P(X).
7) To calculate Q(X), press [D]. The program prompt you for X and then calculate Q(X).
8) To replace P(x) with Q(x), press [E].
9) To calculate a polynomial root, press [f][A]. The program prompts you to enter a guess and a tolerance value. The program assigns a maximum of 50 iterations. The program displays intermediate guess values for the root. If the program finds a real root it displays "RT=<root_value>"

NOTE
I am using B in curly braces instead of square brackets to mean the button B, because this web sites interprets B in square brackets as a directive to display bold text.

Examples

Input of Coefficients for P(X)

Enter the following polynomial (with N=3) by pressing [A]:
P(X) = X^3 - 6*X^2 + 11*X - 6

Values of P(X)

Press [C] to enter any a value for X to calculate P(X). Here is a table of values:

Code:

X    P(X)
--------------
-2    -60
-1    -24
0     -6
0.5  -1.8750
1      0
1.5   0.3750
2      0
2.5  -0.375
3      0
3.5   1.8750
4     6

Values of P'(X)

The first derivative of P(X) is:

P'(X) 3*X^2 - 12*X + 11

Press [f][C] to enter any a value for X to calculate P'(X). Here is a table of values:

Code:

X    P(X)
--------------
-2     47
-1     26
0      11
1        2
2      -1
3       2
4      11

Values of P''(X)

The second derivative of P(X) is:

P''(X) = 6*X - 12

Press [f][D] to enter any a value for X to calculate P''(X). Here is a table of values:

Code:

X    P''(X)
--------------
-2     -24
-1     -18
0      -12
1       -6
2       0
3       6
4       12

Values of Definite Integrals

Press [f][E] to calculate the integral of P(X) between two points. Here is a table of values:

Code:

A    B    Integral
------------------------
0      1      -2.25
0      2      -2
0      3      -2.25
1      2       0.25

Roots of P(X)

Press [f][A] to calculate a real root near a guess you provide. The initial guess of X=4 yields the root of 3. Similarly an initial guess of X=0 yields the root 1.

Deflating and Replacing P(X)

Deflating P(X) with X0=1 yields:

Q(X) = X^2 - 5*X + 6 and R=0, since the value of 1 is a real root of P(X)

Press {B} and enter 1 for the "X0?" prompt. The program displays the above coefficients for Q(X).

Press [E] to replace P(X) with Q(X). You can verify the new coefficients of P(X) by pressing [A} and viewing the new polynomial order and coefficients.

[Thomas Klemm]

Quote:2) To deflate the polynomial P(x), press [B]. The program will prompt you for X0 and then display the coefficients b() of the deflated polynomial Q(x).

It appears that [B] gets interpreted as tag for bold.
But I haven't figured out yet how to correctly escape that.

Cheers
Thomas

[Namir]

I updated the text based on your note, using {B} instead of B in teh square brackets. I also added a note about this issue.

Thanks!