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