# $Source: /u/maple/research/lib/DEtools/diffop/src/RCS/symmetric_power,v $
# $Notify: mvanhoei@daisy.uwaterloo.ca $

# Procedures for users
#    symmetric_power, symmetric_product, exterior_power
# Procedures for use only in the rest of the diffop package:
#    none
# Procedures for use only in this file:
#    symext


########################
#    Various tools     #-----------------------------------------
########################	diffop_tools

macro(
	x=`DEtools/diffop/x`,
	DF=`DEtools/diffop/DF`,

	solve=`solve/linear`,
	symmetric_power=`DEtools/diffop/symmetric_power`,
	symmetric_product=`DEtools/diffop/symmetric_product`,
	symext=`DEtools/diffop/symext`,
	exterior_power=`DEtools/diffop/exterior_power`
):

# Input: f in Qbar(x)[DF] (in RootOf syntax)
# Output: g such that f(y)=0 -> g(y^n)=0
symmetric_power:=proc(f,n)
	global DF,x;
	local i,zero,a,y,ypower,k,subs_rem,co,  L0,L1,A1,A2;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if degree(f,DF)=2 then
		# For order 2 use the method in the paper
		# by Bronstein/Mulders/Weil.
		L0,L1:=1,DF;
		A1:=normal(coeff(f,DF,0)/lcoeff(f,DF));
		A2:=normal(coeff(f,DF,1)/lcoeff(f,DF));
		for i from 1 to n do
			L0,L1:=L1,collect(L1*(DF+i*A2)+diff(L1,x)
				+i*(n-(i-1))*A1*L0,DF,evala@Normal)
		od;
		return L1
	fi;
	ypower:=y(x)^n;
	zero:=a[0]*ypower;
	subs_rem:=diff(y(x),seq(x,k=1..degree(f,DF)))=-convert([y(x)*
	 coeff(f,DF,0)/lcoeff(f,DF),seq(diff(y(x),seq(x,k=1..i))*
	 coeff(f,DF,i)/lcoeff(f,DF),i=1..degree(f,DF)-1)],`+`);
	co:=combinat[numbcomp](degree(f,DF)+n,degree(f,DF));
	symext(co,subs_rem,a,ypower,zero)
end:


# This is not a mathematically correct name, but we use it anyway
# because of the resemblance with symmetric_power.
# Input: a sequence of differential operators (no optional arguments allowed).
symmetric_product:=proc()
	local i,zero,a,ypower,k,subs_rem,co,j,y;
	global DF,x;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	ypower:=convert([seq(y[i](x),i=1..nargs)],`*`);
	zero:=a[0]*ypower;
	subs_rem:={seq(diff(y[j](x),seq(x,k=1..degree(args[j],DF)))=-convert([y[j](x)*
	 coeff(args[j],DF,0)/lcoeff(args[j],DF),seq(diff(y[j](x),seq(x,k=1..i))*
	 coeff(args[j],DF,i)/lcoeff(args[j],DF),i=1..degree(args[j],DF)-1)],`+`)
	,j=1..nargs)};
	co:=convert([seq(degree(i,DF),i=args)],`*`);
	symext(co,subs_rem,a,ypower,zero)
end:

symext:=proc(co,subs_rem,a)
	local i,z,zero,res,k,vars,b;
	global DF,x;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	z:=args[4];
	zero:=args[5];
	vars:={seq(a[k],k=0..co)};
	b[0]:=1;
	for i from 1 to co do
		z:=evala(Normal(subs(subs_rem,diff(z,x))),'independent');
		b[i]:=collect(denom(z)*(diff(b[i-1],x)+b[i-1]*DF),DF,evala@Normal);
		z:=numer(z);
		zero:=expand(zero+a[i]*z)
	od;
	userinfo(5,'diffop',`solving linear equations...`);
	res:=subs(solve({coeffs(zero,{seq(`if`(has(i,indets(args[4],function)),i,NULL)
	 ,i=indets(zero))})},vars),convert([seq(a[k]*b[k],k=0..co)],`+`));
	userinfo(5,'diffop',`done solving linear equations.`);
	z:=nops(indets(res) intersect vars)-1;
	collect(subs(solve({coeff(res,DF,co-z)=1,seq(coeff(res,DF,i)=0
	,i=co-z+1..co)},vars),res),DF,evala@Normal)
end:

# Also called: associated equation
exterior_power:=proc(f,d)
	local a,z,i,j,subs_rem,k,co,zero,y;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if d<=0 then RETURN('undefined') elif d>degree(f,DF) then RETURN(1) fi;
	z:=LinearAlgebra:-LA_Main:-Determinant(Matrix(d,d,[
	[seq(y[i](x),i=1..d)],seq([seq(diff(y[i](x),x$j),i=1..d)],j=1..d-1)]),
         method=LinearAlgebra:-NoUserValue);
	subs_rem:=[seq(
	 diff(y[j](x),seq(x,k=1..degree(f,DF)))=-convert([y[j](x)*
	 coeff(f,DF,0)/lcoeff(f,DF),seq(diff(y[j](x),seq(x,k=1..i))*
	 coeff(f,DF,i)/lcoeff(f,DF),i=1..degree(f,DF)-1)],`+`),j=1..d)];
	co:=binomial(degree(f,DF),d);
	zero:=a[0]*z;
	symext(co,subs_rem,a,z,zero)
end:

#savelib('symmetric_power','symmetric_product','symext','exterior_power'):