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