Polynomial Ops Program - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: HP-41C Software Library (/forum-11.html) +--- Thread: Polynomial Ops Program (/thread-297.html) Polynomial Ops Program - Namir - 01-01-2014 03:01 PM 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*X^N-1) + R Where: b = a b = b * X0 + a 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*X^N-1 P''(X) = 2*a2 + 6*a3*X + 12*a4*X^2 + ... + N*(N-1)*a*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*(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=" 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. RE: Polynomial Ops Program - Thomas Klemm - 01-01-2014 05:50 PM 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 RE: Polynomial Ops Program - Namir - 01-01-2014 07:34 PM 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!