/* * Copyright (c) 1993 David I. Bell * Permission is granted to use, distribute, or modify this source, * provided that this copyright notice remains intact. * * Calculate with polynomials of one variable. */ obj poly {deg, coef}; define pol() { local x, d, i; d = param(0) - 1; if (d < 0) quit "No coefficients for pol"; if (d == 0) return param(1); obj poly x; x.deg = d; mat x.coef[d+1]; for (i = 0; i <= d; i++) x.coef[d-i] = param(i+1); return x; } define poly_print(a) { local i, n; for (i = a.deg; i >= 0; i--) { n = a.coef[i]; if (n == 0) continue; if (i == a.deg) { if (isreal(n) && (n < 0)) { print "- " : ; n = abs(n); } } else { if (!isreal(n) || (n > 0)) print " + " : ; else { print " - " : ; n = abs(n); } } if ((n != 1) && (i > 0)) { if (isreal(n)) print n : "*" : ; else print "(" : n : ")*" : ; } switch (i) { case 0: if (isreal(n)) print n : ; else print "(" : n : ")" : ; break; case 1: print "X" : ; break; default: print "X^" : i : ; } } } define poly_add(a, b) { local x, d; if (isnum(b)) { x = a; x.coef[0] += b; return x; } if (isnum(a)) { x = b; x.coef[0] += a; return x; } if (a.deg == b.deg) { d = a.deg; while (a.coef[d] == -b.coef[d]) if (--d <= 0) return a.coef[0] + b.coef[0]; } d = max(a.deg, b.deg); obj poly x; x.deg = d; mat x.coef[d+1]; while (d >= 0) { if (d > a.deg) x.coef[d] = b.coef[d]; else if (d > b.deg) x.coef[d] = a.coef[d]; else x.coef[d] = a.coef[d] + b.coef[d]; d--; } return x; } define poly_neg(a) { local x, i; x = a; for (i = x.deg; i >= 0; i--) x.coef[i] = -x.coef[i]; return x; } define poly_sub(a, b) { return a + (-b); } define poly_mul(a, b) { local x, i, j; if (isnum(b)) { if (b == 0) return 0; if (b == 1) return a; if (b == -1) return -a; x = a; for (i = x.deg; i >= 0; i--) x.coef[i] *= b; return x; } if (isnum(a)) { if (a == 0) return 0; if (a == 1) return a; if (a == -1) return -a; x = b; for (i = x.deg; i >= 0; i--) x.coef[i] *= a; return x; } obj poly x; x.deg = a.deg + b.deg; mat x.coef[x.deg+1]; for (i = a.deg; i >= 0; i--) for (j = b.deg; j >= 0; j--) x.coef[i+j] += a.coef[i] * b.coef[j]; return x; } define poly_div(a, b) { local i, x; if (!isnum(b)) quit "Only division by numbers currently allowed"; if (b == 0) quit "Division by zero"; if (b == 1) return a; if (b == -1) return -a; x = a; for (i = x.deg; i >= 0; i--) x.coef[i] /= b; return x; } define ev(a, x) { local i, r; obj poly r; if (!istype(a, r)) quit "Evaluating non-polynomial"; i = a.deg; r = a.coef[i]; while (--i >= 0) r = r * x + a.coef[i]; return r; } global lib_debug; if (!isnum(lib_debug) || lib_debug>0) print "obj poly {deg, coef} defined" if (!isnum(lib_debug) || lib_debug>0) print "pol() defined" if (!isnum(lib_debug) || lib_debug>0) print "poly_print(a) defined" if (!isnum(lib_debug) || lib_debug>0) print "poly_add(a, b) defined" if (!isnum(lib_debug) || lib_debug>0) print "poly_neg(a) defined" if (!isnum(lib_debug) || lib_debug>0) print "poly_sub(a, b) defined" if (!isnum(lib_debug) || lib_debug>0) print "poly_mul(a, b) defined" if (!isnum(lib_debug) || lib_debug>0) print "poly_div(a, b) defined" if (!isnum(lib_debug) || lib_debug>0) print "ev(a, x) defined" if (!isnum(lib_debug) || lib_debug>0) print "Use pol() to make polynomials (high coefficient first)" if (!isnum(lib_debug) || lib_debug>0) print "Use ev(a, x) to evaluate them"