##
##                       ####################
##                       #       poly       #
##                       ####################
##
##                     routines for 
##         type checking of polynomials in several variabes,
##         converting polynomials into vectors,
##         solving linear equations between polynomials,
##         determining if a polynomial belongs to the
##             vector space spanned by other polynomials, 
##
##copyleft David Joyner, 6-21-96
#############################################################         

`type/POLYNOMIAL`:=proc(f,v::{name,list(name)})
 type(f,polynom(anything,v));
end:

##### global_n is a global variable representing the
####  number of variables expressing the monomial
`type/MY_VARIABLES`:=proc(f::anything)
global global_n;
member(f,[seq(x.i,i=1..global_n)]);
end:

`type/MONOMIAL`:=proc(f)
 global global_n;
 local good,i;
 good:=`false`;
 if nops(f)=1 and type(f,MY_VARIABLES) then RETURN(`true`); fi;
 if nops(f)>1 then
  for i from 1 to nops(f) do
   if (type(op(i,f),MY_VARIABLES) or (type(op(i,f),`^`) 
       and type(op(1,op(i,f)),MY_VARIABLES))) then
         good:=`true` 
   else
    good:=`false`;
    break;
   fi;
  od; 
 fi;
 RETURN(good);
end:

monomial_part:=proc(f)
##
## f should be a constant times a monomial
##
 global global_n;
 local i,L;
 L:=1;
 for i from 1 to nops(f) do
  if type(op(i,f),MY_VARIABLES) or 
     (type(op(i,f),`^`) and type(op(1,op(i,f)),MY_VARIABLES)) then              
       L:=L*op(i,f); 
  fi;
 od;
 if type(f,`^`) then L:=f; fi; 
 RETURN(L);
end:

##
##the proc below must be modified to deal with polynomials in variables
##other than x1,x2
##
enum_of_monomial_old:=proc(f)
##
##only works for degree 2
##
 global global_n;
local n,d,local_n;
local_n:=2;
if f<>monomial_part(f) then print(`Not a monomial`); RETURN(); fi;
d:=degree(f,[seq(x.i,i=1..local_n)]);
n:=d*(d+1)/2+degree(f,x1)+1;
#### d(d+1)/2 essentially counts the monomials of degree < d in x1, x2
RETURN(n);
end:

partial_binomial_sum:=proc(d::integer)
###
###this counts the monomials of total degree <d in global_n variables 
###
global global_n;
local psum,k;
psum:=sum(binomial(k+global_n-1,global_n-1),k=1..d-1);  ### changed
RETURN(psum);
end:

##
## the following proc returns the number associated
## to a monomial as follows: we order the monomials in x1,...,xn
## first by their total degree then by their order of 
## occurence in (x1+...+xn)^k.
## We call this number the "enumeration" of the monomial
##
##the number of linearly independent terms in (x1+...+xn)^k
## is binomial(k+n-1,n-1)
##
### idea: expand (x1+x2+...+xn)^k (k=1,2,3,...) and 
### enumerate the terms in order of occurence - so the
### n terms for k=1 go first, the terms for k=2 go 2nd, ... 
### for the enumeration of f, read off its total degree
### and then compare f with the monomial part of each term
### 
enum_of_monomial:=proc(f::anything)
global global_n;
local m,d,e,i,a;
if f<>monomial_part(f) then print(`Not a monomial`); RETURN(); fi;
d:=degree(f,[seq(x.i,i=1..global_n)]);
e:=find_in_multi_exp(f);
if d>2 then a:=binomial(d-1,global_n-2); else a:=0; fi; #### changed
m:=partial_binomial_sum(d)+e;
#print(d,e,m);
RETURN(m);
end:

find_in_multi_exp:=proc(f::anything)
global global_n;
local d,i,j,jd1,jd2,f0;
#print(f,` 1 findinmultiexp`);
if f<>monomial_part(f) then print(`Not a monomial`); RETURN(); fi;
d:=degree(f,[seq(x.i,i=1..global_n)]);
#print(d,` 2 findinmultiexp`);
i:='i';
jd1:=sum('x.i',i=1..global_n);
jd2:=expand(jd1^d);
#print(jd1,jd2,` 3 findinmultiexp`);
for i from 1 to nops(jd2) do
  f0[i]:=monomial_part(op(i,jd2));
  if f0[i]=f then 
#   print(f,i);
   RETURN(i); fi;
od;
end:

##
##the proc below is the inverse of the proc above
##
monomial_of_enum:=proc(m::integer)
local p,i,j,d,degree;
degree:=0;
for i from 1 to m do
 if i*(i+1)/2>=m then degree:=i-1; break; fi;
od;
d:=m-degree*(degree+1)/2-1;
#print(` total, x1 degree `,degree,d);
p:=x1^d*x2^(degree-d);
RETURN(p);
end:

##returns list of monomial terms
##
monomial_terms:=proc(f)
local L,i,m;
m:=nops(f);
L:=[];
for i from 1 to m do
 L:=[op(L),monomial_part(op(i,f))];
od;
RETURN(L);
end:

##
##returns list of coefficients
##
coeff_part:=proc(f)
 global global_n;
 local i,m,L,fterms;
m:=nops(f);
L:=[];
for i from 1 to m do
 fterms[i]:=op(i,f);
 if type(op(1,fterms[i]),integer) then
  L:=[op(L),op(1,fterms[i])];
 else
  L:=[op(L),1];
 fi;
od;
RETURN(L);
end:

##
##converts a polynomial in x1,x2 into a vector 
##associated to the monomial terms via the 
##enumeration scheme
##
poly2vector:=proc(f)
local V,v,i,j,m,C,L,max_enum,enum;
m:=nops(f);
C:=coeff_part(f);
L:=monomial_terms(f);
max_enum:=max(seq(enum_of_monomial(L[i]),i=1..m));
V:=n_vectors(max_enum);
v:=[seq(0,i=1..max_enum)];
for i from 1 to m do
 enum:=enum_of_monomial(L[i]);
 v:=v+C[i]*V[enum];
od;
RETURN(v);
end:

n_vectors:=proc(n::integer)
local i,j,e,L;
e[1]:=[1,seq(0,i=1..n-1)];
e[n]:=[seq(0,i=1..n-1),1];
for i from 2 to n-1 do
  e[i]:=[seq(0,j=1..i-1),1,seq(0,j=1..n-i)];
od;
L:=[seq(e[i],i=1..n)];
RETURN(L);
end:

##
##this proc is the inverse to the poly2vector proc
##
vector2poly:=proc(v::list)
local p,i,j,n;
n:=nops(v);
p:=0;
for i from 1 to n do
 p:=p+v[i]*monomial_of_enum(i);
od;
RETURN(p);
end:

##
##this proc will return a basis of the polys
##
poly_basis:=proc(L::list) ##L is a list of polys in x1,x2
local v,vfill,i,j,n,vecs,m,B,b,p,basis_of_polys,enum,max_enum;
n:=nops(L);
max_enum:=max(seq(nops(poly2vector(L[i])),i=1..n));
#print(max_enum);
for i from 1 to n do
 enum[i]:=nops(poly2vector(L[i]));
 vfill[i]:=[op(poly2vector(L[i])),seq(0,j=1..max_enum-enum[i])];
 v[i]:=convert(vfill[i],vector);
od;
vecs:=[seq(v[i],i=1..n)];
B:=linalg[basis](vecs);
m:=nops(B);
for i from 1 to m do
 b[i]:=convert(B[i],list);
 p[i]:=vector2poly(b[i]);
od;
basis_of_polys:=[seq(p[i],i=1..m)];
RETURN(basis_of_polys);
end:

##
##this proc returns 1 if v is in the v space
##spanned by the vectors in B, 0 otherwise
## it is implicitly assumed B is a list of basis vectors
## all vectors are represented as lists
##
is_in:=proc(v::list,B::list)
local i,j,va,vb,A,c,n;
n:=nops(B);
for i from 1 to n do
 vb[i]:=convert(B[i],vector);
od;
va:=convert(v,vector);
A:=linalg[transpose](linalg[matrix]([seq(vb[i],i=1..n)]));
c:=linalg[linsolve](A,va);
if c<>NULL then RETURN(1) else RETURN(0); fi;
end:

is_in_poly:=proc(p,L::list) ##p should be a poly in x1,x2
local i,B,n,max_enum,Ltemp,v,enum,j;
n:=nops(L);
Ltemp[0]:=p;
for i from 1 to n do Ltemp[i]:=L[i]: od;
max_enum:=max(seq(nops(poly2vector(Ltemp[i])),i=0..n));
for i from 0 to n do
 enum[i]:=nops(poly2vector(Ltemp[i]));
 v[i]:=[op(poly2vector(Ltemp[i])),seq(0,j=1..max_enum-enum[i])];
# v[i]:=convert(vfill[i],vector);
od;
#vp:=poly2vector(p);
#for i from 1 to n do
# vB[i]:=poly2vector(L[i]);
#od;
B:=[seq(v[i],i=1..n)];
RETURN(is_in(v[0],B));
end:

##
##the proc tries to find a linear combin of the L[i]'s
##which equals p. Returns 0 if no linear combin exists and
##a list of coeffs otherwise
##
solve_poly:=proc(p,L::list) ##p should be a poly in x1,x2
local i,n,max_enum,Ltemp,v,enum,j,c,va,vb,A,v0;
n:=nops(L);
Ltemp[0]:=p;
for i from 1 to n do Ltemp[i]:=L[i]: od;
max_enum:=max(seq(nops(poly2vector(Ltemp[i])),i=0..n));
for i from 0 to n do
 enum[i]:=nops(poly2vector(Ltemp[i]));
 v[i]:=[op(poly2vector(Ltemp[i])),seq(0,j=1..max_enum-enum[i])];
od;
v0:=v[0];
for i from 1 to n do
 vb[i]:=convert(v[i],vector);
od;
va:=convert(v0,vector);
A:=linalg[transpose](linalg[matrix]([seq(vb[i],i=1..n)]));
c:=linalg[linsolve](A,va);
if c<>NULL then 
 RETURN(convert(c,list)) 
else 
 RETURN(0); 
fi;
end:

#######################################
###### Some differentiation routines
#######################################

mono_diff:=proc(m::anything,g::anything)
##
####   m should be a MONOMIAL (`type/MONOMIAL` may have a bug)
####   g should be a POLYNOMIAL in MY_VARIABLES
##
 global global_n;
 local i,j,k,temp;
 temp:=g;
 k:=nops(m);
 for i from 1 to k do
  if type(op(i,m),`^`) then
   j:=op(2,op(i,m));
  else j:=1;
  fi;
  temp:=diff(temp,x.i$j);
 od;
 RETURN(temp);
end:

###############
# This routine computes the differentiation inner product
# of 2 polynomials in several variables
###############
diff_inner:=proc(f::anything,g::anything)
##
#### f,g should be POLYNOMIALs in MY_VARIABLES
##
 global global_n;
 local j,c,coeffs,m,r,temp,tempf;
 temp:=0;
 if type(f,`+`) then r:=nops(f);
  else r:=0;
 fi;
 if r<>0 then 
  coeffs:=coeff_part(f);
  for j from 1 to r do
   m:=monomial_part(op(j,f));
   c:=coeffs[j];
   temp:=temp+c*mono_diff(m,g);
  od;
 else
  tempf:=monomial_part(f);
  temp:=mono_diff(tempf,g);
 fi;
 RETURN(temp);
end:

save `poly.m`;