poly.cal

Calc Software Package for Number Calc

/*
 * poly - calculate with polynomials of one variable
 *
 * Copyright (C) 1999  Ernest Bowen
 *
 * Calc is open software; you can redistribute it and/or modify it under
 * the terms of the version 2.1 of the GNU Lesser General Public License
 * as published by the Free Software Foundation.
 *
 * Calc is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
 * Public License for more details.
 *
 * A copy of version 2.1 of the GNU Lesser General Public License is
 * distributed with calc under the filename COPYING-LGPL.  You should have
 * received a copy with calc; if not, write to Free Software Foundation, Inc.
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 *
 * @(#) $Revision: 29.2 $
 * @(#) $Id: poly.cal,v 29.2 2000/06/07 14:02:25 chongo Exp $
 * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/poly.cal,v $
 *
 * Under source code control:	1990/02/15 01:50:35
 * File existed as early as:	before 1990
 *
 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
 */

/*
 * A collection of functions designed for calculations involving
 *	polynomials in one variable (by Ernest W. Bowen).
 *
 * On starting the program the independent variable has identifier x
 *	and name "x", i.e. the user can refer to it as x, the
 *	computer displays it as "x".  The name of the independent
 *	variable is stored as varname, so, for example, varname = "alpha"
 *	will change its name to "alpha".  At any time, the independent
 *	variable has only one name.  For some purposes, a name like
 *	"sin(t)" or "(a + b)" or "\lambda" might be useful;
 *	names like "*" or "-27" are legal but might give expressions
 *	that are difficult to intepret.
 *
 * Polynomial expressions may be constructed from numbers and the
 *	independent variable and other polynomials by the algebraic
 *	operations +, -, *, ^, and if the result is a polynomial /.
 *	The operations // and % are defined to have the quotient and
 *	remainder meanings as usually defined for polynomials.
 *
 * When polynomials are assigned to idenfifiers, it is convenient to
 *	think of the polynomials as values.  For example, p = (x - 1)^2
 *	assigns to p a polynomial value in the same way as q = (7 - 1)^2
 *	would assign to q a number value.  As with number expressions
 *	involving operations, the expression used to define the
 *	polynomial is usually lost; in the above example, the normal
 *	computer display for p will be	x^2 - 2x + 1.  Different
 *	identifiers may of course have the same polynomial value.
 *
 * The polynomial we think of as a_0 + a_1 * x + ... + a_n * x^n,
 *	for number coefficients a_0, a_1, ... a_n may also be
 *	constructed as pol(a_0, a_1, ..., a_n).	 Note that here the
 *	coefficients are to be in ascending power order.  The independent
 *	variable is pol(0,1), so to use t, say, as an identifier for
 *	this, one may assign  t = pol(0,1).  To simultaneously specify
 *	an identifier and a name for the independent variable, there is
 *	the instruction var, used as in identifier = var(name).	 For
 *	example, to use "t" in the way "x" is initially, one may give
 *	the instruction	 t = var("t").
 *
 * There are four parameters pmode, order, iod and ims for controlling
 *	the format in which polynomials are displayed.
 *	The parameter pmode may have values "alg" or "list": the
 *	former gives a display as an algebraic formula, while
 *	the latter only lists the coefficients.	 Whether the terms or
 *	coefficients are in ascending or descending power order is
 *	controlled by order being "up" or "down".  If the
 *	parameter iod (for integer-only display), the polynomial
 *	is expressed in terms of a polynomial whose coefficients are
 *	integers with gcd = 1, the leading coefficient having positive
 *	real part, with where necessary a leading multiplying integer,
 *	a Gaussian integer multiplier if the coefficients are complex
 *	with a common complex factor, and a trailing divisor integer.
 *	If a non-zero value is assigned to the parameter ims,
 *	multiplication signs will be inserted where appropriate;
 *	this may be useful if the expression is to be copied to a
 *	program or a string to be used with eval.
 *
 * For evaluation of polynomials the standard function is ev(p, t).
 *	If p is a polynomial and t anything for which the relevant
 *	operations can be performed, this returns the value of p
 *	at t.  The function ev(p, t) also accepts lists or matrices
 *	as possible values for p; each element of p is then evaluated
 *	at t.  For other p, t is ignored and the value of p is returned.
 *	If an identifier, a, say, is used for the polynomial, list or
 *	matrix p, the definition
 *			define a(t) = ev(a, t);
 *	permits a(t) to be used for the value of a at t as if the
 *	polynomial, list or matrix were a function.  For example,
 *	if a = 1 + x^2, a(2) will return the value 5, just as if
 *			define a(t) = 1 + t^2;
 *	had been used.	 However, when the polynomial definition is
 *	used, changing the polynomial a will change a(t) to the value
 *	of the new polynomial at t.  For example,
 *	after
 *		L = list(x, x^2, x^3, x^4);
		define a(t) = ev(a, t);
 *	the loop
 *		for (i = 0; i < 4; i++)
 *			print ev(L[[i]], 5);
 *	may be replaced by
 *		for (i = 0; i < 4; i++) {
 *			a = L[[i]];
 *			print a(5);
 *		}
 *
 * Matrices with polynomial elements may be added, subtracted and
 *	multiplied as long as the usual rules for compatibility are
 *	observed.  Also, matrices may be multiplied by polynomials,
 *	i.e. if p is a	polynomial and A a matrix whose elements
 *	may be numbers or polynomials, p * A returns the matrix of
 *	the same shape as A with each element multiplied by p.
 *	Square matrices may also be 'substituted for the variable' in
 *	polynomials, e.g. if A is an m x m matrix, and
 *	p = x^2 + 3 * x + 2, ev(p, A) returns the same as
 *	A^2 + 3 * A + 2 * I, where I is the unit m x m matrix.
 *
 * On starting this program, three demonstration polynomials a, b, c
 *	have been defined.  The functions a(t), b(t), c(t) corresponding
 *	to a, b, c, and x(t) corresponding to x, have also been
 *	defined, so the usual function notation can be used for
 *	evaluations of a, b, c and x.  For x, as long as x identifies
 *	the independent variable, x(t) should return the value of t,
 *	i.e. it acts as an identity function.
 *
 * Functions defined include:
 *
 *	monic(a) returns the monic multiple of a, i.e., if a != 0,
 *		the multiple of a with leading coefficient 1
 *	conj(a) returns the complex conjugate of a
 *	ispmult(a,b) returns 1 or 0 according as a is or is not
 *		a polynomial multiple of b
 *	pgcd(a,b) returns the monic gcd of a and b
 *	pfgcd(a,b) returns a list of three polynomials (g, u, v)
 *		where g = pgcd(a,b) and g = u * a + v * b.
 *	plcm(a,b) returns the monic lcm of a and b
 *
 *	interp(X,Y,t) returns the value at t of the polynomial given
 *		by Newtonian divided difference interpolation, where
 *		X is a list of x-values, Y a list of corresponding
 *		y-values.  If t is omitted, the interpolating
 *		polynomial is returned.	 A y-value may be replaced by
 *		list (y, y_1, y_2, ...), where y_1, y_2, ... are
 *		the reduced derivatives at the corresponding x;
 *		i.e. y_r is the r-th derivative divided by fact(r).
 *	mdet(A) returns the determinant of the square matrix A,
 *		computed by an algorithm that does not require
 *		inverses;  the built-in det function usually fails
 *		for matrices with polynomial elements.
 *	D(a,n) returns the n-th derivative of a; if n is omitted,
 *		the first derivative is returned.
 *
 * A first-time user can see what the initially defined polynomials
 *	a, b and c are, and experiment with the algebraic operations
 *	and other functions that have been defined by giving
 *	instructions like:
 *			a
 *			b
 *			c
 *			(x^2 + 1) * a
 *			a^27
 *			a * b
 *			a % b
 *			a // b
 *			a(1 + x)
 *			a(b)
 *			conj(c)
 *			g = pgcd(a, b)
 *			g
 *			a / g
 *			D(a)
 *			mat A[2,2] = {1 + x, x^2, 3, 4*x}
 *			mdet(A)
 *			D(A)
 *			A^2
 *			define A(t) = ev(A, t)
 *			A(2)
 *			A(1 + x)
 *			define L(t) = ev(L, t)
 *			L = list(x, x^2, x^3, x^4)
 *			L(5)
 *			a(L)
 *			interp(list(0,1,2,3), list(2,3,5,7))
 *			interp(list(0,1,2), list(0,list(1,0),2))
 *
 * One check on some of the functions is provided by the Cayley-Hamilton
 *	theorem:  if A is any m x m matrix and I the m x m unit matrix,
 *	and x is pol(0,1),
 *			ev(mdet(x * I - A), A)
 *	should return the zero m x m matrix.
 */


obj poly {p};

define pol() {
	local u,i,s;
	obj poly u;
	s = list();
	for (i=1; i<= param(0); i++) append (s,param(i));
	i=size(s) -1;
	while (i>=0 && s[[i]]==0) {i--; remove(s)}
	u.p = s;
	return u;
}

define ispoly(a) {
	local y;
	obj poly y;
	return istype(a,y);
}

define findlist(a) {
	if (ispoly(a)) return a.p;
	if (a) return list(a);
	return list();
}

pmode = "alg";	/* The other acceptable pmode is "list" */
ims = 0;	/* To be non-zero if multiplication signs to be inserted */
iod = 0;	/* To be non-zero for integer-only display */
order = "down"	/* Determines order in which coefficients displayed */

define poly_print(a) {
	local f, g, t;
	if (size(a.p) == 0) {
		print 0:;
		return;
	}
	if (iod) {
		g = gcdcoeffs(a);
		t = a.p[[size(a.p) - 1]] / g;
		if (re(t) < 0) { t = -t; g = -g;}
		if (g != 1) {
			if (!isreal(t)) {
				if (im(t) > re(t)) g *= 1i;
				else if (im(t) <= -re(t)) g *= -1i;
			}
			if (isreal(g)) f = g;
			else f = gcd(re(g), im(g));
			if (num(f) != 1) {
				print num(f):;
				if (ims) print"*":;
			}
			if (!isreal(g)) {
				printf("(%d)", g/f);
				if (ims) print"*":;
			}
			if (pmode == "alg") print"(":;
			polyprint(1/g * a);
			if (pmode == "alg") print")":;
			if (den(f) > 1) print "/":den(f):;
			return;
		}
	}
	polyprint(a);
}

define polyprint(a) {
	local s,n,i,c;
	s = a.p;
	n=size(s) - 1;
	if (pmode=="alg") {
		if (order == "up") {
			i = 0;
			while (!s[[i]]) i++;
			pterm (s[[i]], i);
			for (i++ ; i <= n; i++) {
				c = s[[i]];
				if (c) {
					if (isreal(c)) {
						if (c > 0) print" + ":;
						else {
							print" - ":;
							c = -c;
						}
					}
					else print " + ":;
					pterm(c,i);
				}
			}
			return;
		}
		if (order == "down") {
			pterm(s[[n]],n);
			for (i=n-1; i>=0; i--) {
				c = s[[i]];
				if (c) {
					if (isreal(c)) {
						if (c > 0) print" + ":;
						else {
							print" - ":;
							c = -c;
						}
					}
					else print " + ":;
					pterm(c,i);
				}
			}
			return;
		}
		quit "order to be up or down";
	}
	if (pmode=="list") {
		plist(s);
		return;
	}
	print pmode,:"is unknown mode";
}


define poly_neg(a) {
	local s,i,y;
	obj poly y;
	s = a.p;
	for (i=0; i< size(s); i++) s[[i]] = -s[[i]];
	y.p = s;
	return y;
}

define poly_conj(a) {
	local s,i,y;
	obj poly y;
	s = a.p;
	for (i=0; i < size(s); i++) s[[i]] = conj(s[[i]]);
	y.p = s;
	return y;
}

define poly_inv(a) = pol(1)/a;	/* This exists only for a of zero degree */

define poly_add(a,b) {
	local sa, sb, i, y;
	obj poly y;
	sa=findlist(a); sb=findlist(b);
	if (size(sa) > size(sb)) swap(sa,sb);
	for (i=0; i< size(sa); i++) sa[[i]] += sb[[i]];
	while (i < size(sb)) append (sa, sb[[i++]]);
	while (i > 0 && sa[[--i]]==0) remove (sa);
	y.p = sa;
	return y;
}

define poly_sub(a,b) {
	 return a + (-b);
}

define poly_cmp(a,b) {
	local sa, sb;
	sa = findlist(a);
	sb=findlist(b);
	return	(sa != sb);
}