4.4BSD/usr/src/contrib/calc-1.26.4/lib/poly.cal

/*
 * 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"