poly.cal

Calc Software Package for Number Calc

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 d, u, i, m, n, sa, sb, sq;
	local obj poly q;
	local obj poly 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);

if (config("resource_debug") & 3) {
	print "obj poly {p} defined";
}