#$Id: pde.1.txt,v 2.16 1995/10/09 19:23:03 coffey Exp $
#
#	homoMap
#
# Constructs all monomials of degree g:  z^a w^b Z^c W^d, a+b+c+d = g 
# input	degree g
# output set of coefficient lists [ [1,,,,0]], [1,,,,0], ]
homoMap := proc(g)   
local a,b,c;
{seq(seq(seq([1,a,b,c,g-a-b-c,0],c=0..g-a-b),b=0..g-a),a=0..g)};
end:

#
#	h
#
# Applies group action to a specific coefficient list.
# input: [1,,,,0], p and q picked such that (p,q)=1 and q!=+-1(mod p)
#        typical examples are 2,1 5,2 5,3   Good one: 8,3
# output: a number, if an integer then monomial is invariant under
#         group action.
h := proc(l,p,q)
(-l[2]-l[3]*q+l[4]+l[5]*q)/p;
end:

#
#	Invmonos
#
# Applies the group action:, h, to each item in set and removes those that
# are not invariant under the action.
# input: { [,,,], [,,,], }, p, q   takes set input.
# output: { [,,,], [,,,], }
Invmonos := proc(L,p,q)
local l,x;
l := NULL;
for x in L do
	if type(h(x,p,q),integer) then l := l,x; fi;
od;
{l};
end:

#
#	listdiff
#
# This will take the derivative of every elment in the list with respect
# to the list element specified.
# input: [ [e,a,b,c,d,u], ], x an integer
# output: [ [e,a,b,c,d,u], ]
listdiff := proc(L,x)
local i,l;
	l := NULL;
	for i in L do
		if i[x] > 0 then
			l := l,[i[1]*i[x],op(2..x-1,i),i[x]-1,op(x+1..nops(i),i)];
		fi;
	od;
	[l];
end:
#
#	listprod
#
# This will take two polynomials in list form and take their product.
# no collection is done for I don't know what I will need to collect
# on.
# input: [ [,,,,,], [,,,,,], ], [ [,,,,,], [,,,,,], ]
# output: [ [,,,,,], [,,,,,], ]
# This multiplies the first element, adds the next four, and puts the
# last element in a list.
listprod := proc(L1,L2)
local i,j;
[seq(seq([
		L1[i][1]*L2[j][1],
		L1[i][2]+L2[j][2],
		L1[i][3]+L2[j][3],
		L1[i][4]+L2[j][4],
		L1[i][5]+L2[j][5],
		sort([L1[i][6],L2[j][6]])
		],j=1..nops(L2)),i=1..nops(L1))];
end:


#
#       Rightside
#
# Outputs the right hand side of the quation in list form.
# input: equation number n, degree g
# output: [ [,,,], [,,,], ]
Rightside := proc(n,g)
local k;
if n = 1 then
[seq([1,g-2-k,k,g-k,k,0],k=0..g-2)];
elif n = 2 then
[seq([1,g-2-k,k,g-2-k,k+2,0],k=0..g-2)];
elif n = 3 then
[seq([1,g-1-k,k,g-1-k,k,0],k=0..g-1)];
elif n = 4 then
[seq([1,g-1-k,k,g-1-k,k,0],k=0..g-1)];
elif n = 5 then
[seq([1,g-2-k,k,g-1-k,k+1,0],k=0..g-2)];
elif n = 6 then
[seq([1,g-k-1,k,g-2-k,k+1,0],k=0..g-2)];
fi;
end:

#
#	CoeffGrab
#
# Note: the input to this should be collected on the monomials first.
#
# This procedure will take in two polynomials and a choice monomial
# Each monomial in the first polynomial is then taken and a conjugate
# monomial is constructed that would add to the choice monomial.  This
# monomial is then checked for in the other side.  This is continued
# until one is found in the second polynomial.  It is then multiplied
# and added to a queue.  This queue is then output as a polynomial
#
# This will now also take in the Rightside monomials and subtract them 
# at the end from the current coefficients, if necessary.

CoeffGrab := proc(A,B,M,R)
local m,n,N,G,L,queue,this,j;
	L := NULL;
	for m in A do
# printf(`m = %a\n`,m);
		N := [1,M[2]-m[2],M[3]-m[3],M[4]-m[4],M[5]-m[5],0];
		G := NULL;
		for this in B do
			if [op(2..5,this)] = [op(2..5,N)] then G := G,this; fi;	
		od;
		G := [G];
		queue := NULL;
		for n in G do
			queue := queue,[n[1]*m[1],sort([n[6],m[6]])];
		od;	
		if queue <> NULL then L := L,queue; fi;
	od;
	for j from 1 to nops(R) do
		if op(2..5,M) = op(2..5,R[j]) then 
			L := L,[-R[j][1],[0,0]];
		fi;
	od;
	[L];
end:

#
#	RightPoly
#
# This constructs the righthand side in list form.
# It also changes the sign on the constant term so that they
# can just be appended to the list of monomials.
#
RightPoly := proc(n,g)
local k;

if n = 1 then
 [seq([(3*g/(4*g+8)-1/4)*(g-2)!/((g-k-2)!*k!),g-k-2,k,g-k,k,0],k=0..g-2)];
elif n = 2 then
 [seq([(3*g/(4*g+8)-1/4)*(g-2)!/((g-k-2)!*k!),g-k-2,k,g-k-2,k+2,0],k=0..g-2)];
elif n = 3 then
 [[1/2,0,g-1,0,g-1,0],seq([((g-2)!/((g-2-k)!*k!))*((g-1)/(2*(g-1-k))+3*g/(4*g+8)-1/4),g-1-k,k,g-1-k,k,0],k=0..g-2)];
elif n = 4 then
 [seq([1/2*(g-2)!/((g-2-k)!*k!),g-1-k,k,g-1-k,k,0],k=0..g-2),
  seq([(3*g/(4*g+8)+1/4)*(g-2)!/((g-2-k)!*k!),g-2-k,k+1,g-2-k,k+1,0],k=0..g-2)];
elif n = 5 then
 [seq([(3*g/(4*g+8)-1/4)*(g-2)!/((g-2-k)!*k!),g-2-k,k,g-1-k,k+1,0],k=0..g-2)];
elif n = 6 then
 [seq([(3*g/(4*g+8)-1/4)*(g-2)!/((g-2-k)!*k!),g-1-k,k,g-2-k,k+1,0],k=0..g-2)];
fi;
end:

#
#       RightSide
#
# This proc outpus the polynomials on the right hand side.
# input: equation number n, degree g
# output: [ [,,,], [,,,], ]
RightSide := proc(n,g)
local k;
if n = 1 then
[seq([1,g-2-k,k,g-k,k,0],k=0..g-2)];
elif n = 2 then
[seq([1,g-2-k,k,g-2-k,k+2,0],k=0..g-2)];
elif n = 3 then
[seq([1,g-1-k,k,g-1-k,k,0],k=0..g-1)];
elif n = 4 then
[seq([1,g-1-k,k,g-1-k,k,0],k=0..g-1)];
elif n = 5 then
[seq([1,g-2-k,k,g-1-k,k+1,0],k=0..g-2)];
elif n = 6 then
[seq([1,g-k-1,k,g-2-k,k+1,0],k=0..g-2)];
fi;
end:


#
#	CollectMonomials
#
# This proc should take in a polynomial in list form and collect on all
# the monomials by adding together the constant coefficients.  This will
# only occur if a,b,c,d, and the u term are the same.
#
#CollectMonomials := proc(A);
#local m,B;
#	B := NULL;
#	while nops(A) > 0 do 
#		m := A[1];
#		A := [op(2..nops(A),A)];
#		for i from 1 to nops(A) do
#			if [op(2..6,A[i])] = [op(2..6,m)] then
#				 m := [m[1]+A[i][1],op(2..6,m)];
#				 
#			fi;
#		od;
#	od;
#
#end:
# **********************************************************
#
#	basic
#
basic := proc(g,p,q)
local homo,n,harmonics,nakedMono,polyhold,a,b,c,d,capN,i,coeff1,coeff2,
	  mono1,mono2,f,fz,fw,fZ,fW,arbitrary,fzfz,fzfZ,fwfw,fwfW,fzfw,fwfZ,
      cl,RHS,st,T;

homo := homoMap(g);
n := Invmonos(homo,p,q);

print(`Constructing harmonics`);
harmonics := NULL;
arbitrary := 0;
while nops(n) > 0 do
#print(nops(n));	
    arbitrary := arbitrary+1;
	nakedMono := n[1];
#print(`NakedMono`,nakedMono);
	n := n minus { nakedMono };
	polyhold := [op(1..5,nakedMono),arbitrary];
	a := nakedMono[2];	b := nakedMono[3];
	c := nakedMono[4];	d := nakedMono[5];

	capN := max(min(a,c),min(b,d));
	for i from 0 to capN do
#print(capN,i);
		coeff1:=(-1)^(i+1)*product((a-k)*(c-k)/((b+k+1)*(d+k+1)),k=0..i);		
#print(`Coeff1`,coeff1);
		mono1:=[coeff1,a-i-1,b+i+1,c-i-1,d+i+1,arbitrary];
		coeff2:=(-1)^(i+1)*product((b-k)*(d-k)/((a+k+1)*(c+k+1)),k=0..i);
#print(`Coeff2`,coeff2);
		mono2:=[coeff2,a+i+1,b-i-1,c+i+1,d-i-1,arbitrary];
#print(`Mono1: `,mono1,` Mono2: `,mono2);
	
		if coeff1 <> 0 then polyhold := polyhold,mono1; fi;
		if coeff2 <> 0 then polyhold := polyhold,mono2; fi;
		n := n minus { [1,a-i-1,b+i+1,c-i-1,d+i+1,0], [1,a+i+1,b-i-1,c+i+1,d-i-1,0] };
	od;
#print(`Polyhold`,polyhold);
	harmonics := harmonics,polyhold;
od;
f := [harmonics];

print(`Starting differentiation`);
fz := listdiff(f,2); 
fw := listdiff(f,3);
fZ := listdiff(f,4);
fW := listdiff(f,5);

save(f,fz,fw,fZ,fW,`data_diff.txt`);

print(`Starting product and colleciton.`);
cl := NULL;
print(`Polynomial #1`);
RHS := Rightside(1,g);
st := time();
cl := cl,seq(CoeffGrab(fz,fz,RHS[i],RightPoly(1,g)),i=1..nops(RHS));
print(time()-st);
T := `Poly 1 done!`; save(T,`cl.txt`);
print(`Polynomial #2`);
RHS := Rightside(2,g);
st := time();
cl := cl,seq(CoeffGrab(fw,fw,RHS[i],RightPoly(2,g)),i=1..nops(RHS));
print(time()-st);
T := `Poly 2 done!`; save(T,`cl.txt`);
print(`Polynomial #3`);
RHS := Rightside(3,g);
st := time();
cl := cl,seq(CoeffGrab(fz,fZ,RHS[i],RightPoly(3,g)),i=1..nops(RHS));
print(time()-st);
T := `Poly 3 done!`; save(T,`cl.txt`);
print(`Polynomial #4`);
RHS := Rightside(4,g);
st := time();
cl := cl,seq(CoeffGrab(fw,fW,RHS[i],RightPoly(4,g)),i=1..nops(RHS));
print(time()-st);
T := `Poly 4 done!`; save(T,`cl.txt`);
print(`Polynomial #5`);
RHS := Rightside(5,g);
st := time();
cl := cl,seq(CoeffGrab(fz,fw,RHS[i],RightPoly(5,g)),i=1..nops(RHS));
print(time()-st);
T := `Poly 5 done!`; save(T,`cl.txt`);
print(`Polynomial #6`);
RHS := Rightside(6,g);
st := time();
cl := cl,seq(CoeffGrab(fw,fZ,RHS[i],RightPoly(6,g)),i=1..nops(RHS));
print(time()-st);
T := `Poly 6 done!`; save(T,`cl.txt`);

#equ := {op(EndListToPoly([cl]))};
#solve(equ,indets(equ,function));

print(`Finished constructing equations.`);
{cl};

#[f,fz,fw,fZ,fW];

# Now all we need to do is solve the system generated by cl.  That is,
# we need to convert the cl notation to something solvable by, e.g. matlab.


#print(`Starting products`);
#fzfz := listprod([fz],[fz]);
#fzfZ := listprod([fz],[fZ]);
#fwfw := listprod([fw],[fw]);
#fwfW := listprod([fw],[fW]);
#fzfw := listprod([fz],[fw]);
#fwfZ := listprod([fw],[fZ]);
#
#save(fzfz,fzfZ,fwfw,fwfW,fzfw,fwfZ,`data_prod.txt`);

end:


EndListToPoly := proc(L)
local i,j;
subs(u_0_0=1,{seq(convert([seq(L[i][j][1]*cat(u,_,convert(L[i][j][2][1],string),_,convert(L[i][j][2][2],string)),j=1..nops(L[i]))],`+`),i=1..nops(L))});
end:

PolyListToPoly := proc(L)
local i,j;
convert([seq(cat(u,_,convert(L[i][6],string))*L[i][1]*z^L[i][2]*w^L[i][3]*zbar^L[i][4]*wbar^L[i][5],i=1..nops(L))],`+`);
end: