1/* 2 * poly - calculate with polynomials of one variable 3 * 4 * Copyright (C) 1999,2021 Ernest Bowen 5 * 6 * Calc is open software; you can redistribute it and/or modify it under 7 * the terms of the version 2.1 of the GNU Lesser General Public License 8 * as published by the Free Software Foundation. 9 * 10 * Calc is distributed in the hope that it will be useful, but WITHOUT 11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 12 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General 13 * Public License for more details. 14 * 15 * A copy of version 2.1 of the GNU Lesser General Public License is 16 * distributed with calc under the filename COPYING-LGPL. You should have 17 * received a copy with calc; if not, write to Free Software Foundation, Inc. 18 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 19 * 20 * Under source code control: 1990/02/15 01:50:35 21 * File existed as early as: before 1990 22 * 23 * Share and enjoy! :-) http://www.isthe.com/chongo/tech/comp/calc/ 24 */ 25 26/* 27 * A collection of functions designed for calculations involving 28 * polynomials in one variable (by Ernest W. Bowen). 29 * 30 * On starting the program the independent variable has identifier x 31 * and name "x", i.e. the user can refer to it as x, the 32 * computer displays it as "x". The name of the independent 33 * variable is stored as varname, so, for example, varname = "alpha" 34 * will change its name to "alpha". At any time, the independent 35 * variable has only one name. For some purposes, a name like 36 * "sin(t)" or "(a + b)" or "\lambda" might be useful; 37 * names like "*" or "-27" are legal but might give expressions 38 * that are difficult to interpret. 39 * 40 * Polynomial expressions may be constructed from numbers and the 41 * independent variable and other polynomials by the algebraic 42 * operations +, -, *, ^, and if the result is a polynomial /. 43 * The operations // and % are defined to have the quotient and 44 * remainder meanings as usually defined for polynomials. 45 * 46 * When polynomials are assigned to identifiers, it is convenient to 47 * think of the polynomials as values. For example, p = (x - 1)^2 48 * assigns to p a polynomial value in the same way as q = (7 - 1)^2 49 * would assign to q a number value. As with number expressions 50 * involving operations, the expression used to define the 51 * polynomial is usually lost; in the above example, the normal 52 * computer display for p will be x^2 - 2x + 1. Different 53 * identifiers may of course have the same polynomial value. 54 * 55 * The polynomial we think of as a_0 + a_1 * x + ... + a_n * x^n, 56 * for number coefficients a_0, a_1, ... a_n may also be 57 * constructed as pol(a_0, a_1, ..., a_n). Note that here the 58 * coefficients are to be in ascending power order. The independent 59 * variable is pol(0,1), so to use t, say, as an identifier for 60 * this, one may assign t = pol(0,1). To simultaneously specify 61 * an identifier and a name for the independent variable, there is 62 * the instruction var, used as in identifier = var(name). For 63 * example, to use "t" in the way "x" is initially, one may give 64 * the instruction t = var("t"). 65 * 66 * There are four parameters pmode, order, iod and ims for controlling 67 * the format in which polynomials are displayed. 68 * The parameter pmode may have values "alg" or "list": the 69 * former gives a display as an algebraic formula, while 70 * the latter only lists the coefficients. Whether the terms or 71 * coefficients are in ascending or descending power order is 72 * controlled by order being "up" or "down". If the 73 * parameter iod (for integer-only display), the polynomial 74 * is expressed in terms of a polynomial whose coefficients are 75 * integers with gcd = 1, the leading coefficient having positive 76 * real part, with where necessary a leading multiplying integer, 77 * a Gaussian integer multiplier if the coefficients are complex 78 * with a common complex factor, and a trailing divisor integer. 79 * If a non-zero value is assigned to the parameter ims, 80 * multiplication signs will be inserted where appropriate; 81 * this may be useful if the expression is to be copied to a 82 * program or a string to be used with eval. 83 * 84 * For evaluation of polynomials the standard function is ev(p, t). 85 * If p is a polynomial and t anything for which the relevant 86 * operations can be performed, this returns the value of p 87 * at t. The function ev(p, t) also accepts lists or matrices 88 * as possible values for p; each element of p is then evaluated 89 * at t. For other p, t is ignored and the value of p is returned. 90 * If an identifier, a, say, is used for the polynomial, list or 91 * matrix p, the definition 92 * define a(t) = ev(a, t); 93 * permits a(t) to be used for the value of a at t as if the 94 * polynomial, list or matrix were a function. For example, 95 * if a = 1 + x^2, a(2) will return the value 5, just as if 96 * define a(t) = 1 + t^2; 97 * had been used. However, when the polynomial definition is 98 * used, changing the polynomial a will change a(t) to the value 99 * of the new polynomial at t. For example, 100 * after 101 * L = list(x, x^2, x^3, x^4); 102 define a(t) = ev(a, t); 103 * the loop 104 * for (i = 0; i < 4; i++) 105 * print ev(L[[i]], 5); 106 * may be replaced by 107 * for (i = 0; i < 4; i++) { 108 * a = L[[i]]; 109 * print a(5); 110 * } 111 * 112 * Matrices with polynomial elements may be added, subtracted and 113 * multiplied as long as the usual rules for compatibility are 114 * observed. Also, matrices may be multiplied by polynomials, 115 * i.e. if p is a polynomial and A a matrix whose elements 116 * may be numbers or polynomials, p * A returns the matrix of 117 * the same shape as A with each element multiplied by p. 118 * Square matrices may also be 'substituted for the variable' in 119 * polynomials, e.g. if A is an m x m matrix, and 120 * p = x^2 + 3 * x + 2, ev(p, A) returns the same as 121 * A^2 + 3 * A + 2 * I, where I is the unit m x m matrix. 122 * 123 * On starting this program, three demonstration polynomials a, b, c 124 * have been defined. The functions a(t), b(t), c(t) corresponding 125 * to a, b, c, and x(t) corresponding to x, have also been 126 * defined, so the usual function notation can be used for 127 * evaluations of a, b, c and x. For x, as long as x identifies 128 * the independent variable, x(t) should return the value of t, 129 * i.e. it acts as an identity function. 130 * 131 * Functions defined include: 132 * 133 * monic(a) returns the monic multiple of a, i.e., if a != 0, 134 * the multiple of a with leading coefficient 1 135 * conj(a) returns the complex conjugate of a 136 * ispmult(a,b) returns 1 or 0 according as a is or is not 137 * a polynomial multiple of b 138 * pgcd(a,b) returns the monic gcd of a and b 139 * pfgcd(a,b) returns a list of three polynomials (g, u, v) 140 * where g = pgcd(a,b) and g = u * a + v * b. 141 * plcm(a,b) returns the monic lcm of a and b 142 * 143 * interp(X,Y,t) returns the value at t of the polynomial given 144 * by Newtonian divided difference interpolation, where 145 * X is a list of x-values, Y a list of corresponding 146 * y-values. If t is omitted, the interpolating 147 * polynomial is returned. A y-value may be replaced by 148 * list (y, y_1, y_2, ...), where y_1, y_2, ... are 149 * the reduced derivatives at the corresponding x; 150 * i.e. y_r is the r-th derivative divided by fact(r). 151 * mdet(A) returns the determinant of the square matrix A, 152 * computed by an algorithm that does not require 153 * inverses; the built-in det function usually fails 154 * for matrices with polynomial elements. 155 * D(a,n) returns the n-th derivative of a; if n is omitted, 156 * the first derivative is returned. 157 * 158 * A first-time user can see what the initially defined polynomials 159 * a, b and c are, and experiment with the algebraic operations 160 * and other functions that have been defined by giving 161 * instructions like: 162 * a 163 * b 164 * c 165 * (x^2 + 1) * a 166 * a^27 167 * a * b 168 * a % b 169 * a // b 170 * a(1 + x) 171 * a(b) 172 * conj(c) 173 * g = pgcd(a, b) 174 * g 175 * a / g 176 * D(a) 177 * mat A[2,2] = {1 + x, x^2, 3, 4*x} 178 * mdet(A) 179 * D(A) 180 * A^2 181 * define A(t) = ev(A, t) 182 * A(2) 183 * A(1 + x) 184 * define L(t) = ev(L, t) 185 * L = list(x, x^2, x^3, x^4) 186 * L(5) 187 * a(L) 188 * interp(list(0,1,2,3), list(2,3,5,7)) 189 * interp(list(0,1,2), list(0,list(1,0),2)) 190 * 191 * One check on some of the functions is provided by the Cayley-Hamilton 192 * theorem: if A is any m x m matrix and I the m x m unit matrix, 193 * and x is pol(0,1), 194 * ev(mdet(x * I - A), A) 195 * should return the zero m x m matrix. 196 */ 197 198 199obj poly {p}; 200 201define pol() { 202 local u,i,s; 203 obj poly u; 204 s = list(); 205 for (i=1; i<= param(0); i++) append (s,param(i)); 206 i=size(s) -1; 207 while (i>=0 && s[[i]]==0) {i--; remove(s)} 208 u.p = s; 209 return u; 210} 211 212define ispoly(a) { 213 local y; 214 obj poly y; 215 return istype(a,y); 216} 217 218define findlist(a) { 219 if (ispoly(a)) return a.p; 220 if (a) return list(a); 221 return list(); 222} 223 224pmode = "alg"; /* The other acceptable pmode is "list" */ 225ims = 0; /* To be non-zero if multiplication signs to be inserted */ 226iod = 0; /* To be non-zero for integer-only display */ 227order = "down" /* Determines order in which coefficients displayed */ 228 229define poly_print(a) { 230 local f, g, t; 231 if (size(a.p) == 0) { 232 print 0:; 233 return; 234 } 235 if (iod) { 236 g = gcdcoeffs(a); 237 t = a.p[[size(a.p) - 1]] / g; 238 if (re(t) < 0) { t = -t; g = -g;} 239 if (g != 1) { 240 if (!isreal(t)) { 241 if (im(t) > re(t)) g *= 1i; 242 else if (im(t) <= -re(t)) g *= -1i; 243 } 244 if (isreal(g)) f = g; 245 else f = gcd(re(g), im(g)); 246 if (num(f) != 1) { 247 print num(f):; 248 if (ims) print"*":; 249 } 250 if (!isreal(g)) { 251 printf("(%d)", g/f); 252 if (ims) print"*":; 253 } 254 if (pmode == "alg") print"(":; 255 polyprint(1/g * a); 256 if (pmode == "alg") print")":; 257 if (den(f) > 1) print "/":den(f):; 258 return; 259 } 260 } 261 polyprint(a); 262} 263 264define polyprint(a) { 265 local s,n,i,c; 266 s = a.p; 267 n=size(s) - 1; 268 if (pmode=="alg") { 269 if (order == "up") { 270 i = 0; 271 while (!s[[i]]) i++; 272 pterm (s[[i]], i); 273 for (i++ ; i <= n; i++) { 274 c = s[[i]]; 275 if (c) { 276 if (isreal(c)) { 277 if (c > 0) print" + ":; 278 else { 279 print" - ":; 280 c = -c; 281 } 282 } 283 else print " + ":; 284 pterm(c,i); 285 } 286 } 287 return; 288 } 289 if (order == "down") { 290 pterm(s[[n]],n); 291 for (i=n-1; i>=0; i--) { 292 c = s[[i]]; 293 if (c) { 294 if (isreal(c)) { 295 if (c > 0) print" + ":; 296 else { 297 print" - ":; 298 c = -c; 299 } 300 } 301 else print " + ":; 302 pterm(c,i); 303 } 304 } 305 return; 306 } 307 quit "order to be up or down"; 308 } 309 if (pmode=="list") { 310 plist(s); 311 return; 312 } 313 print pmode,:"is unknown mode"; 314} 315 316 317define poly_neg(a) { 318 local s,i,y; 319 obj poly y; 320 s = a.p; 321 for (i=0; i< size(s); i++) s[[i]] = -s[[i]]; 322 y.p = s; 323 return y; 324} 325 326define poly_conj(a) { 327 local s,i,y; 328 obj poly y; 329 s = a.p; 330 for (i=0; i < size(s); i++) s[[i]] = conj(s[[i]]); 331 y.p = s; 332 return y; 333} 334 335define poly_inv(a) = pol(1)/a; /* This exists only for a of zero degree */ 336 337define poly_add(a,b) { 338 local sa, sb, i, y; 339 obj poly y; 340 sa=findlist(a); sb=findlist(b); 341 if (size(sa) > size(sb)) swap(sa,sb); 342 for (i=0; i< size(sa); i++) sa[[i]] += sb[[i]]; 343 while (i < size(sb)) append (sa, sb[[i++]]); 344 while (i > 0 && sa[[--i]]==0) remove (sa); 345 y.p = sa; 346 return y; 347} 348 349define poly_sub(a,b) { 350 return a + (-b); 351} 352 353define poly_cmp(a,b) { 354 local sa, sb; 355 sa = findlist(a); 356 sb=findlist(b); 357 return (sa != sb); 358} 359 360define poly_mul(a,b) { 361 local sa,sb,i, j, y; 362 if (ismat(a)) swap(a,b); 363 if (ismat(b)) { 364 y = b; 365 for (i=matmin(b,1); i <= matmax(b,1); i++) 366 for (j = matmin(b,2); j<= matmax(b,2); j++) 367 y[i,j] = a * b[i,j]; 368 return y; 369 } 370 obj poly y; 371 sa=findlist(a); sb=findlist(b); 372 y.p = listmul(sa,sb); 373 return y; 374} 375 376define listmul(a,b) { 377 local da,db, s, i, j, u; 378 da=size(a)-1; db=size(b)-1; 379 s=list(); 380 if (da >= 0 && db >= 0) { 381 for (i=0; i<= da+db; i++) { u=0; 382 for (j = max(0,i-db); j <= min(i, da); j++) 383 u += a[[j]]*b[[i-j]]; append (s,u);}} 384 return s; 385} 386 387define ev(a,t) { 388 local v, i, j; 389 if (ismat(a)) { 390 v = a; 391 for (i = matmin(a,1); i <= matmax(a,1); i++) 392 for (j = matmin(a,2); j <= matmax(a,2); j++) 393 v[i,j] = ev(a[i,j], t); 394 return v; 395 } 396 if (islist(a)) { 397 v = list(); 398 for (i = 0; i < size(a); i++) 399 append(v, ev(a[[i]], t)); 400 return v; 401 } 402 if (!ispoly(a)) return a; 403 if (islist(t)) { 404 v = list(); 405 for (i = 0; i < size(t); i++) 406 append(v, ev(a, t[[i]])); 407 return v; 408 } 409 if (ismat(t)) return evpm(a.p, t); 410 return evp(a.p, t); 411} 412 413define evp(s,t) { 414 local n,v,i; 415 n = size(s); 416 if (!n) return 0; 417 v = s[[n-1]]; 418 for (i = n - 2; i >= 0; i--) v=t * v +s[[i]]; 419 return v; 420} 421 422define evpm(s,t) { 423 local m, n, V, i, I; 424 n = size(s); 425 m = matmax(t,1) - matmin(t,1); 426 if (matmax(t,2) - matmin(t,2) != m) quit "Non-square matrix"; 427 mat V[m+1, m+1]; 428 if (!n) return V; 429 mat I[m+1, m+1]; 430 matfill(I, 0, 1); 431 V = s[[n-1]] * I; 432 for (i = n - 2; i >= 0; i--) V = t * V + s[[i]] * I; 433 return V; 434} 435pzero = pol(0); 436x = pol(0,1); 437varname = "x"; 438define x(t) = ev(x, t); 439 440define iszero(a) { 441 if (ispoly(a)) 442 return !size(a.p); 443 return a == 0; 444} 445 446define isstring(a) = istype(a, " "); 447 448define var(name) { 449 if (!isstring(name)) quit "Argument of var is to be a string"; 450 varname = name; 451 return pol(0,1); 452} 453 454define pcoeff(a) { 455 if (isreal(a)) print a:; 456 else print "(":a:")":; 457} 458 459define pterm(a,n) { 460 if (n==0) { 461 pcoeff(a); 462 return; 463 } 464 if (n==1) { 465 if (a!=1) { 466 pcoeff(a); 467 if (ims) print"*":; 468 } 469 print varname:; 470 return; 471 } 472 if (a!=1) { 473 pcoeff(a); 474 if (ims) print"*":; 475 } 476 print varname:"^":n:; 477} 478 479define plist(s) { 480 local i, n; 481 n = size(s); 482 print "( ":; 483 if (order == "up") { 484 for (i=0; i< n-1 ; i++) 485 print s[[i]]:",",:; 486 if (n) print s[[i]],")":; 487 else print "0 )":; 488 } 489 else { 490 if (n) print s[[n-1]]:; 491 for (i = n - 2; i >= 0; i--) 492 print ", ":s[[i]]:; 493 print " )":; 494 } 495} 496 497define deg(a) = size(a.p) - 1; 498 499define polydiv(a,b) { 500 local d, u, i, m, n, sa, sb, sq; 501 local obj poly q; 502 local obj poly r; 503 sa=findlist(a); sb = findlist(b); sq = list(); 504 m=size(sa)-1; n=size(sb)-1; 505 if (n<0) quit "Zero divisor"; 506 if (m<n) return list(pzero, a); 507 d = sb[[n]]; 508 while ( m >= n) { u = sa[[m]]/d; 509 for (i = 0; i< n; i++) sa[[m-n+i]] -= u*sb[[i]]; 510 push(sq,u); remove(sa); m--; 511 while (m>=n && sa[[m]]==0) { m--; remove(sa); push(sq,0)}} 512 while (m>=0 && sa[[m]]==0) { m--; remove(sa);} 513 q.p = sq; r.p = sa; 514 return list(q, r);} 515 516define poly_mod(a,b) { 517 local u; 518 u=polydiv(a,b); 519 return u[[1]]; 520} 521 522define poly_quo(a,b) { 523 local p; 524 p = polydiv(a,b); 525 return p[[0]]; 526} 527 528define ispmult(a,b) = iszero(a % b); 529 530define poly_div(a,b) { 531 if (!ispmult(a,b)) quit "Result not a polynomial"; 532 return poly_quo(a,b); 533} 534 535define pgcd(a,b) { 536 local r; 537 if (iszero(a) && iszero(b)) return pzero; 538 while (!iszero(b)) { 539 r = a % b; 540 a = b; 541 b = r; 542 } 543 return monic(a); 544} 545 546define plcm(a,b) = monic( a * b // pgcd(a,b)); 547 548define pfgcd(a,b) { 549 local u, v, u1, v1, s, q, r, d, w; 550 u = v1 = pol(1); v = u1 = pol(0); 551 while (size(b.p) > 0) {s = polydiv(a,b); 552 q = s[[0]]; 553 a = b; b = s[[1]]; u -= q*u1; v -= -q*v1; 554 swap(u,u1); swap(v,v1);} 555 d=size(a.p)-1; if (d>=0 && (w= 1/a.p[[d]]) !=1) 556 { a *= w; u *= w; v *= w;} 557 return list(a,u,v); 558} 559 560define monic(a) { 561 local s, c, i, d, y; 562 if (iszero(a)) return pzero; 563 obj poly y; 564 s = findlist(a); 565 d = size(s)-1; 566 for (i=0; i<=d; i++) s[[i]] /= s[[d]]; 567 y.p = s; 568 return y; 569} 570 571define coefficient(a,n) = (n < size(a.p)) ? a.p[[n]] : 0; 572 573define D(a, n) { 574 local i,j,v; 575 if (isnull(n)) n = 1; 576 if (!isint(n) || n < 1) quit "Bad order for derivative"; 577 if (ismat(a)) { 578 v = a; 579 for (i = matmin(a,1); i <= matmax(a,1); i++) 580 for (j = matmin(a,2); j <= matmax(a,2); j++) 581 v[i,j] = D(a[i,j], n); 582 return v; 583 } 584 if (!ispoly(a)) return 0; 585 return Dp(a,n); 586} 587 588define Dp(a,n) { 589 local i, v; 590 if (n > 1) return Dp(Dp(a, n-1), 1); 591 obj poly v; 592 v.p=list(); 593 for (i=1; i<size(a.p); i++) append (v.p, i*a.p[[i]]); 594 return v; 595} 596 597 598define cgcd(a,b) { 599 if (isreal(a) && isreal(b)) return gcd(a,b); 600 while (a) { 601 b -= bround(b/a) * a; 602 swap(a,b); 603 } 604 if (re(b) < 0) b = -b; 605 if (im(b) > re(b)) b *= -1i; 606 else if (im(b) <= -re(b)) b *= 1i; 607 return b; 608} 609 610define gcdcoeffs(a) { 611 local s,i,g, c; 612 s = a.p; 613 g=0; 614 for (i=0; i < size(s) && g != 1; i++) 615 if (c = s[[i]]) g = cgcd(g, c); 616 return g; 617} 618 619define interp(X, Y, t) = evalfd(makediffs(X,Y), t); 620 621define makediffs(X,Y) { 622 local U, D, d, x, y, i, j, k, m, n, s; 623 U = D = list(); 624 n = size(X); 625 if (size(Y) != n) quit"Arguments to be lists of same size"; 626 for (i = n-1; i >= 0; i--) { 627 x = X[[i]]; 628 y = Y[[i]]; 629 m = size(U); 630 if (isnum(y)) { 631 d = y; 632 for (j = 0; j < m; j++) { 633 d = D[[j]] = (D[[j]]-d)/(U[[j]] - x); 634 } 635 push(U, x); 636 push(D, y); 637 } 638 else { 639 s = size(y); 640 for (k = 0; k < s ; k++) { 641 d = y[[k]]; 642 for (j = 0; j < m; j++) { 643 d = D[[j]] = (D[[j]] - d)/(U[[j]] - x); 644 } 645 } 646 for (j=s-1; j >=0; j--) { 647 push(U,x); 648 push(D, y[[j]]); 649 } 650 } 651 } 652 return list(U, D); 653} 654 655define evalfd(T, t) { 656 local U, D, n, i, v; 657 if (isnull(t)) t = pol(0,1); 658 U = T[[0]]; 659 D = T[[1]]; 660 n = size(U); 661 v = D[[n-1]]; 662 for (i = n-2; i >= 0; i--) 663 v = v * (t - U[[i]]) + D[[i]]; 664 return v; 665} 666 667 668define mdet(A) { 669 local n, i, j, k, I, J; 670 n = matmax(A,1) - (i = matmin(A,1)); 671 if (matmax(A,2) - (j = matmin(A,2)) != n) 672 quit "Non-square matrix for mdet"; 673 I = J = list(); 674 k = n + 1; 675 while (k--) { 676 append(I,i++); 677 append(J,j++); 678 } 679 return M(A, n+1, I, J); 680} 681 682define M(A, n, I, J) { 683 local v, J0, i, j, j1; 684 if (n == 1) return A[ I[[0]], J[[0]] ]; 685 v = 0; 686 i = remove(I); 687 for (j = 0; j < n; j++) { 688 J0 = J; 689 j1 = delete(J0, j); 690 v += (-1)^(n-1+j) * A[i, j1] * M(A, n-1, I, J0); 691 } 692 return v; 693} 694 695define mprint(A) { 696 local i,j; 697 if (!ismat(A)) quit "Argument to be a matrix"; 698 for (i = matmin(A,1); i <= matmax(A,1); i++) { 699 for (j = matmin(A,2); j <= matmax(A,2); j++) 700 printf("%8.4d ", A[i,j]); 701 printf("\n"); 702 } 703} 704 705obj poly a; 706obj poly b; 707obj poly c; 708 709define a(t) = ev(a,t); 710define b(t) = ev(b,t); 711define c(t) = ev(c,t); 712 713a=pol(1,4,4,2,3,1); 714b=pol(5,16,8,1); 715c=pol(1+2i,3+4i,5+6i); 716 717if (config("resource_debug") & 3) { 718 print "obj poly {p} defined"; 719} 720