projects / unix-history / blob
? search:
summary | tags | clone url | shortlog | log | commit | commitdiff | tree
blame | history | raw | HEAD
BSD 4_4_Lite2 development
[unix-history] / usr / src / contrib / calc-2.9.3t6 / lib / poly.cal
/*
* 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);
}
define poly_mul(a,b) {
local sa,sb,i, j, y;
if (ismat(a)) swap(a,b);
if (ismat(b)) {
y = b;
for (i=matmin(b,1); i <= matmax(b,1); i++)
for (j = matmin(b,2); j<= matmax(b,2); j++)
y[i,j] = a * b[i,j];
return y;
}
obj poly y;
sa=findlist(a); sb=findlist(b);
y.p = listmul(sa,sb);
return y;
}
define listmul(a,b) {
local da,db, s, i, j, u;
da=size(a)-1; db=size(b)-1;
s=list();
if (da >= 0 && db >= 0) {
for (i=0; i<= da+db; i++) { u=0;
for (j = max(0,i-db); j <= min(i, da); j++)
u += a[[j]]*b[[i-j]]; append (s,u);}}
return s;
}
define ev(a,t) {
local v, i, j;
if (ismat(a)) {
v = a;
for (i = matmin(a,1); i <= matmax(a,1); i++)
for (j = matmin(a,2); j <= matmax(a,2); j++)
v[i,j] = ev(a[i,j], t);
return v;
}
if (islist(a)) {
v = list();
for (i = 0; i < size(a); i++)
append(v, ev(a[[i]], t));
return v;
}
if (!ispoly(a)) return a;
if (islist(t)) {
v = list();
for (i = 0; i < size(t); i++)
append(v, ev(a, t[[i]]));
return v;
}
if (ismat(t)) return evpm(a.p, t);
return evp(a.p, t);
}
define evp(s,t) {
local n,v,i;
n = size(s);
if (!n) return 0;
v = s[[n-1]];
for (i = n - 2; i >= 0; i--) v=t * v +s[[i]];
return v;
}
define evpm(s,t) {
local m, n, V, i, I;
n = size(s);
m = matmax(t,1) - matmin(t,1);
if (matmax(t,2) - matmin(t,2) != m) quit "Non-square matrix";
mat V[m+1, m+1];
if (!n) return V;
mat I[m+1, m+1];
matfill(I, 0, 1);
V = s[[n-1]] * I;
for (i = n - 2; i >= 0; i--) V = t * V + s[[i]] * I;
return V;
}
pzero = pol(0);
x = pol(0,1);
varname = "x";
define x(t) = ev(x, t);
define iszero(a) {
if (ispoly(a))
return !size(a.p);
return a == 0;
}
define isstring(a) = istype(a, " ");
define var(name) {
if (!isstring(name)) quit "Argument of var is to be a string";
varname = name;
return pol(0,1);
}
define pcoeff(a) {
if (isreal(a)) print a:;
else print "(":a:")":;
}
define pterm(a,n) {
if (n==0) {
pcoeff(a);
return;
}
if (n==1) {
if (a!=1) {
pcoeff(a);
if (ims) print"*":;
}
print varname:;
return;
}
if (a!=1) {
pcoeff(a);
if (ims) print"*":;
}
print varname:"^":n:;
}
define plist(s) {
local i, n;
n = size(s);
print "( ":;
if (order == "up") {
for (i=0; i< n-1 ; i++)
print s[[i]]:",",:;
if (n) print s[[i]],")":;
else print "0 )":;
}
else {
if (n) print s[[n-1]]:;
for (i = n - 2; i >= 0; i--)
print ", ":s[[i]]:;
print " )":;
}
}
define deg(a) = size(a.p) - 1;
define polydiv(a,b) {
local q, r, d, u, i, m, n, sa, sb, sq;
obj poly q, r;
sa=findlist(a); sb = findlist(b); sq = list();
m=size(sa)-1; n=size(sb)-1;
if (n<0) quit "Zero divisor";
if (m<n) return list(pzero, a);
d = sb[[n]];
while ( m >= n) { u = sa[[m]]/d;
for (i = 0; i< n; i++) sa[[m-n+i]] -= u*sb[[i]];
push(sq,u); remove(sa); m--;
while (m>=n && sa[[m]]==0) { m--; remove(sa); push(sq,0)}}
while (m>=0 && sa[[m]]==0) { m--; remove(sa);}
q.p = sq; r.p = sa;
return list(q, r);}
define poly_mod(a,b) {
local u;
u=polydiv(a,b);
return u[[1]];
}
define poly_quo(a,b) {
local p;
p = polydiv(a,b);
return p[[0]];
}
define ispmult(a,b) = iszero(a % b);
define poly_div(a,b) {
if (!ispmult(a,b)) quit "Result not a polynomial";
return poly_quo(a,b);
}
define pgcd(a,b) {
local r;
if (iszero(a) && iszero(b)) return pzero;
while (!iszero(b)) {
r = a % b;
a = b;
b = r;
}
return monic(a);
}
define plcm(a,b) = monic( a * b // pgcd(a,b));
define pfgcd(a,b) {
local u, v, u1, v1, s, q, r, d, w;
u = v1 = pol(1); v = u1 = pol(0);
while (size(b.p) > 0) {s = polydiv(a,b);
q = s[[0]];
a = b; b = s[[1]]; u -= q*u1; v -= -q*v1;
swap(u,u1); swap(v,v1);}
d=size(a.p)-1; if (d>=0 && (w= 1/a.p[[d]]) !=1)
{ a *= w; u *= w; v *= w;}
return list(a,u,v);
}
define monic(a) {
local s, c, i, d, y;
if (iszero(a)) return pzero;
obj poly y;
s = findlist(a);
d = size(s)-1;
for (i=0; i<=d; i++) s[[i]] /= s[[d]];
y.p = s;
return y;
}
define coefficient(a,n) = (n < size(a.p)) ? a.p[[n]] : 0;
define D(a, n) {
local i,j,v;
if (isnull(n)) n = 1;
if (!isint(n) || n < 1) quit "Bad order for derivative";
if (ismat(a)) {
v = a;
for (i = matmin(a,1); i <= matmax(a,1); i++)
for (j = matmin(a,2); j <= matmax(a,2); j++)
v[i,j] = D(a[i,j], n);
return v;
}
if (!ispoly(a)) return 0;
return Dp(a,n);
}
define Dp(a,n) {
local i, v;
if (n > 1) return Dp(Dp(a, n-1), 1);
obj poly v;
v.p=list();
for (i=1; i<size(a.p); i++) append (v.p, i*a.p[[i]]);
return v;
}
define cgcd(a,b) {
if (isreal(a) && isreal(b)) return gcd(a,b);
while (a) {
b -= bround(b/a) * a;
swap(a,b);
}
if (re(b) < 0) b = -b;
if (im(b) > re(b)) b *= -1i;
else if (im(b) <= -re(b)) b *= 1i;
return b;
}
define gcdcoeffs(a) {
local s,i,g, c;
s = a.p;
g=0;
for (i=0; i < size(s) && g != 1; i++)
if (c = s[[i]]) g = cgcd(g, c);
return g;
}
define interp(X, Y, t) = evalfd(makediffs(X,Y), t);
define makediffs(X,Y) {
local U, D, d, x, y, i, j, k, m, n, s;
U = D = list();
n = size(X);
if (size(Y) != n) quit"Arguments to be lists of same size";
for (i = n-1; i >= 0; i--) {
x = X[[i]];
y = Y[[i]];
m = size(U);
if (isnum(y)) {
d = y;
for (j = 0; j < m; j++) {
d = D[[j]] = (D[[j]]-d)/(U[[j]] - x);
}
push(U, x);
push(D, y);
}
else {
s = size(y);
for (k = 0; k < s ; k++) {
d = y[[k]];
for (j = 0; j < m; j++) {
d = D[[j]] = (D[[j]] - d)/(U[[j]] - x);
}
}
for (j=s-1; j >=0; j--) {
push(U,x);
push(D, y[[j]]);
}
}
}
return list(U, D);
}
define evalfd(T, t) {
local U, D, n, i, v;
if (isnull(t)) t = pol(0,1);
U = T[[0]];
D = T[[1]];
n = size(U);
v = D[[n-1]];
for (i = n-2; i >= 0; i--)
v = v * (t - U[[i]]) + D[[i]];
return v;
}
define mdet(A) {
local n, i, j, k, I, J;
n = matmax(A,1) - (i = matmin(A,1));
if (matmax(A,2) - (j = matmin(A,2)) != n)
quit "Non-square matrix for mdet";
I = J = list();
k = n + 1;
while (k--) {
append(I,i++);
append(J,j++);
}
return M(A, n+1, I, J);
}
define M(A, n, I, J) {
local v, J0, i, j, j1;
if (n == 1) return A[ I[[0]], J[[0]] ];
v = 0;
i = remove(I);
for (j = 0; j < n; j++) {
J0 = J;
j1 = delete(J0, j);
v += (-1)^(n-1+j) * A[i, j1] * M(A, n-1, I, J0);
}
return v;
}
define mprint(A) {
local i,j;
if (!ismat(A)) quit "Argument to be a matrix";
for (i = matmin(A,1); i <= matmax(A,1); i++) {
for (j = matmin(A,2); j <= matmax(A,2); j++)
printf("%8.4d ", A[i,j]);
printf("\n");
}
}
obj poly a;
obj poly b;
obj poly c;
define a(t) = ev(a,t);
define b(t) = ev(b,t);
define c(t) = ev(c,t);
a=pol(1,4,4,2,3,1);
b=pol(5,16,8,1);
c=pol(1+2i,3+4i,5+6i);
global lib_debug;
if (lib_debug >= 0) {
print "obj poly {p} defined";
print "pol() defined";
print "poly_print(a) defined";
print "poly_add(a, b) defined";
print "poly_sub(a, b) defined";
print "poly_mul(a, b) defined";
print "poly_div(a, b) defined";
print "poly_quo(a,b) defined";
print "poly_mod(a,b) defined";
print "poly_neg(a) defined";
print "poly_conj(a) defined";
print "poly_cmp(a,b) defined";
print "iszero(a) defined";
print "plist(a) defined";
print "listmul(a,b) defined";
print "ev(a,t) defined";
print "evp(s,t) defined";
print "ispoly(a) defined";
print "isstring(a) defined";
print "var(name) defined";
print "pcoeff(a) defined";
print "pterm(a,n) defined";
print "deg(a) defined";
print "polydiv(a,b) defined";
print "D(a,n) defined";
print "Dp(a,n) defined";
print "pgcd(a,b) defined";
print "plcm(a,b) defined";
print "monic(a) defined";
print "pfgcd(a,b) defined";
print "interp(X,Y,x) defined";
print "makediffs(X,Y) defined";
print "evalfd(T,x) defined";
print "mdet(A) defined";
print "M(A,n,I,J) defined";
print "mprint(A) defined";
}
Continuous source code history from original PDP-7 UNIX to FreeBSD 2, the first fully unencumbered FreeBSD release.