#################################################################
#    Package to aid in making generalized qDyson conjectures    #
#           Andrew V. Sills, Georgia Southern Univ              #
#################################################################

PressOn:=false:
Verbose:=false:
Debug:=false:

print(`Generalized qDyson conjecture package`);
print(`by A.V. Sills, Georgia Southern University`);
print(`Version of May 22, 2007`);
print(`accompanies the paper "Disturbing the q-Dyson Conjecture"`);
print(`submitted to the Conference Proceedings`);
print(`for the special session in Experimental Math`);
print(`held at the 2007 AMS-MAA Joint Meetings in New Orleans`);
print(` `);
print(`Enter 'ez()' for a list of procedures`);

ez:=proc();
print(` `);
print(`qSF(a,q,n)....q-shifted factorial`);
print(`qDysonProd(x,a)`);
print(`qMultinomial(a1,a2,...,an)`);
print(`esp1(a,n)...first elementary symmetric polynomial`);
print(`c1m1(r,s, [a1,a2,...,an])`);
print(`c2m1m1(r,s,t,[a1,a2,...,an])`);
print(`c11m1(r,s,t,u,[a1,a2,...,an])`);
print(`Conj1m1(r,s, n)`);
print(`Conj2m1m1(r,s,t, n)`);
print(`Conj11m1m1(r,s,t,u, n)`);
print(`Global boolean variables: Verbose, Debug`):
end:

qSF:=proc(a,q,n) local i; product(1-a*q^i,i=0..n-1) end:

qDysonProd:=proc(x,a) local i,j,n;
  n:=nops(a);
  product(product(
    qSF(x[i]*q/x[j],q,a[j]) *qSF(x[j]/x[i],q,a[i]), j=i+1..n),i=1..n-1)
end:

qMultinomial:=proc() local i,n,a;
  a:=[args];
  n:=nargs;
  qSF(q,q, add(a[i],i=1..n))/product( qSF(q,q,a[i]),i=1..n)
end:

esp1:=proc(a,n) local i; add(a[i],i=1..n) end:

CT:=proc(p,x,n) local i,tmp;
tmp:=p;
for i from 1 to n do
  tmp:= coeff(tmp,x[i],0)
od;
tmp
end:

c1m1:=proc(r,s, a) local n;
n:=nops(a);
CT( 
expand(
  normal(
    coeff(coeff(expand(qDysonProd(x,a)),x[r],1),x[s],-1)
    /(-qMultinomial(op(a)) * (1-q^a[s])  / (1-q^(1+esp1(a,n)-a[s]))
    )
  )
)
,x,n)
end:

c2m1m1:=proc(r,s,t,a) local n;
n:=nops(a);
CT(
expand(
  normal( 
    coeff(coeff(coeff(expand(qDysonProd(x,a)),x[r],2),x[s],-1),x[t],-1)
    /(qMultinomial(op(a)) * (1-q^a[s]) * (1-q^a[t]) / (1-q^(1+esp1(a,n)-a[s]))/
     (1-q^(1+esp1(a,n)-a[t]))/(1-q^(1+esp1(a,n)-a[s]-a[t])))
  )
)
,x,n) 
end:

c11m1m1:=proc(r,s,t,u, a) local n;
n:=nops(a);
CT(
expand(
  normal(
    coeff(coeff(coeff(coeff(
        expand(qDysonProd(x,a)),x[r],1),x[s],1),x[t],-1),x[u],-1)
    /(qMultinomial(op(a)) * (1-q^a[t]) * (1-q^a[u]) / (1-q^(1+esp1(a,n)-a[t]))/
     (1-q^(1+esp1(a,n)-a[u]))/(1-q^(1+esp1(a,n)-a[u]-a[t])))
  )
)
,x,n)
end:


LMdecomp:=proc(p,F1,F2) local p1;
if type(p,`+`) and nops(p)=2 then      
  #return([SKIP,SKIP])
  p1:=sort(p); 
  if Debug then print([degree(op(2,p1),q),degree(op(1,p1),q)-degree(F2,q)]) fi:
  return([ degree(op(2,p1),q), degree(op(1,p1),q) - degree(F2,q)]) 
fi;
if not type(p,`+`) or nops(p)<>4
then print(`Input must be a four term polynomial in q; instead got`,p) 
fi;
p1:=sort(p);
if expand(op(-1,p1)*F1 + op(-2,p1)*F2)=p 
   then return([degree(op(-1,p1),q), degree(op(-2,p1),q)])
elif expand(op(-2,p1)*F1 + op(-1,p1)*F2)=p
   then return([degree(op(-2,p1),q), degree(op(-1,p1),q)])
else ERROR(`The decompostion has failed.`)
fi
end:

DataVectors:=proc(n,restr) local i,k,offset,tmp,zerov; 
zerov:=[seq(0,i=1..n)]; 
tmp:=[seq( [seq(0,i=1..k-1), 1, seq(0,i=k+1..n)], k=1..n),zerov];
offset:=[seq(0,i=1..n)]+ 
      add( [seq(0,i=1..op(k,restr)-1),1,seq(0,i=op(k,restr)+1..n)], 
        k=1..nops(restr));
[seq( tmp[i]+offset, i=1..nops(tmp))]
end:

Conj1m1:=proc(r,s, n)  local a,i,j,k,L,M,S,x,Leq,Lvar,newLHS,
  tmpSubs,newLeq,Lsoln,LRHS ;
a:=[a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12];
a:=[seq(a[i],i=1..n)];
S:=DataVectors(n,[s]);
Lvar:={seq(lambda[i],i=0..n)};
Leq:={};
L:= lambda[0] + add(lambda[k]*a[k],k=1..n);
for j from 1 to n+1 do
  tmpSubs:={seq(a[k]=S[j][k],k=1..n)};
  if Debug then  print(tmpSubs) fi;
  if Debug then  print(c1m1(r,s, S[j] )) fi;
  LRHS:= degree( c1m1(r,s, S[j] ), q); 
  newLeq := subs(tmpSubs, L)=LRHS;
  if Debug then  print(newLeq) fi;
  Leq:= Leq union {newLeq};
od;
Lsoln:= solve(Leq,Lvar); 
if Verbose then
  print(`The coefficient of`, x[r]/x[s],` in the n=`,n);
  print(`case of the q-Dyson product is conjectured to be`);
  print(
  -q^(subs(Lsoln,L))* (1-q^a[s]) / (1-q^(1+esp1(a,n)-a[s])) *
``(q,q)[esp1(a,n)]/product(``(q,q)[a[i]],i=1..n) 
 )
else print([r,s],subs(Lsoln,L))
fi
end:

Check1m1:=proc(n) local c,data,i,j;
c:=combinat[permute](n,2):
data:=combinat[permute](10,n):
for i from 1 to nops(c) do
  print(`------------------------------------------`);
  print( [r,s] =c[i] );   
  for j from 1 to nops(data) do
    if Check1m1One( op(c[i]), data[j])=1 then
      if Verbose then print( data[j],`ok`) fi
    else ERROR(`Problem encountered at [r,s]`=[r,s],`...data value`,data[j])
    fi
  od;
  print(`OK `)       
od
end: 

Check1m1One:=proc(r,s,a) local k,L,n;
n:=nops(a);
if   r<s then L:=1+esp1(a,n)-add(a[k],k=r..s) 
elif r>s then L:=add(a[k],k=s+1..r-1)
else ERROR(`r and s must be distinct integers`)
fi:
#if a[s]=0 then return(0) else
 simplify(
  CT(
   expand(
    normal(
      coeff(
       coeff(
        expand(
          qDysonProd(x,a)
        ),
       x[r],1),
      x[s],-1) 
       *(1-q^(1+esp1(a,n)-a[s]))
       /(-q^L * qMultinomial(op(a)) * (1-q^a[s])  )
      )
     )
  ,x,n)
  )
#fi
end:

Conj2m1m1:=proc(r,s,t,n) local a,i,j,k,L,M,S,x,Leq,Meq,Lvar,Mvar,newLHS,
  CL,CMd,tmpSubs,bothRHS,newLeq, newMeq, Lsoln,Msoln,LRHS,MRHS;
a:=[a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12];
a:=[seq(a[i],i=1..n)];
S:=DataVectors(n,[s,t]);
Lvar:={seq(lambda[i],i=0..n)};
Mvar:={seq(mu[i],i=0..n)};
Leq:={}; Meq:={};
for j from 1 to n+1 do
  tmpSubs:={seq(a[k]=S[j][k],k=1..n)};
#  print(tmpSubs);
#  print(c2m1m1(r,s,t, subs(tmpSubs,a)));
  bothRHS:= LMdecomp( c2m1m1(r,s,t, subs(tmpSubs,a)) ,
             subs(tmpSubs, 1-q^(1+esp1(a,n))), 
             subs(tmpSubs,1-q^(1+esp1(a,n)-a[s]-a[t])) );
  LRHS:=bothRHS[1]; MRHS:=bothRHS[2]; 
  newLeq := subs(tmpSubs, lambda[0] + add(a[k]*lambda[k], k=1..n))=LRHS;
  newMeq := subs(tmpSubs, mu[0] + add(a[k]*mu[k],k=1..n))=MRHS;
#  print(newLeq,newMeq);
  Leq:= Leq union {newLeq};
  Meq:= Meq union {newMeq}
od;
Lsoln:= solve(Leq,Lvar);  
#print(Lsoln);
Msoln:= solve(Meq,Mvar);
L:= lambda[0] + add(lambda[k]*a[k],k=1..n);
M:= mu[0] + add(mu[k]*a[k],k=1..n);
#print(Msoln);
if r<s and s<t 
    then CL:=2+2*esp1(a,n)-2*add(a[k],k=r..t)+add(a[k],k=s+1..t-1);
         CMd:= a[t]
elif s<r and r<t 
    #then CL:= 1  + esp1(a,n)-add(a[k],k=r..t)+ 2*add(a[k],k=s+1..r-1) 
    #     - add(a[k],k=s..r-1);
    then CL:= 1  + esp1(a,n)-add(a[k],k=s..t)+ 2*add(a[k],k=s+1..r-1);
         CMd:= a[s]
elif s<t and t<r 
    then CL:= 2*add(a[k],k=t+1..r-1) + add(a[k],k=s+1..t-1);
         CMd:= a[t]
fi;

if Verbose then
  print(`The coefficient of`, x[r]^2/x[s]/x[t],` in the n=`,n);
  print(`case of the q-Dyson product`);
  #print(qDysonProd(x,a));
  print(` is conjectured to be`);
  print(
  ``(q,q)[esp1(a,n)]/product(``(q,q)[a[i]],i=1..n) * (1-q^(a[s])) * 
   (1-q^(a[t]))
   /(1-q^(1+esp1(a,n)-a[s]))/(1-q^(1+esp1(a,n)-a[t]))/(1-q^(1+esp1(a,n)-a[s]-a[t]))
   *( q^(subs(Lsoln,L)) * (1-q^(1+esp1(a,n) )) 
   + q^(subs(Msoln,M)) * (1-q^(1+esp1(a,n)-a[s]-a[t]))    )
   )
# else print([r,s,t],[subs(Lsoln,L), subs(Msoln,M)],subs(Msoln,M)-subs(Lsoln,L) )
else print([r,s,t],[subs(Lsoln,L), subs(Msoln,M)-subs(Lsoln,L)],`...`,
            [CL-subs(Lsoln,L), CMd-subs(Msoln,M)+subs(Lsoln,L)]
                       )


fi
end:
 
Conj11m1m1:=proc(r,s,t,u,  n) local a,i,j,k,L,M,S,x,Leq,Meq,Lvar,Mvar,newLHS,
  tmpSubs,bothRHS,newLeq, newMeq, Lsoln,Msoln,LRHS,MRHS,CL,CMd;
a:=[a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12];
a:=[seq(a[i],i=1..n)];
S:=DataVectors(n,[t,u]);
Lvar:={seq(lambda[i],i=0..n)};
Mvar:={seq(mu[i],i=0..n)};
Leq:={}; Meq:={};
for j from 1 to n+1 do
  tmpSubs:={seq(a[k]=S[j][k],k=1..n)};
  if Debug then  print(tmpSubs) fi;
  if Debug then  print(c11m1m1(r,s,t,u, subs(tmpSubs,a)),`...`, 
          q^L * (1-q^(1+subs(tmpSubs,esp1(a,n))))
         +q^M * (1-q^(1+subs(tmpSubs,esp1(a,n)-a[t]-a[u])))) 
  fi;
  bothRHS:= LMdecomp( c11m1m1(r,s,t,u, subs(tmpSubs,a)) ,
             subs(tmpSubs, 1-q^(1+esp1(a,n))),
             subs(tmpSubs,1-q^(1+esp1(a,n)-a[t]-a[u])) );
  if bothRHS=[SKIP,SKIP] then 
     print([r,s,t,u], `Result is inconclusive, best to skip this case`); 
     if not PressOn then return() fi 
  fi;
  LRHS:=bothRHS[1]; MRHS:=bothRHS[2];
  newLeq := subs(tmpSubs, lambda[0] + add(a[k]*lambda[k], k=1..n))=LRHS;
  newMeq := subs(tmpSubs, mu[0] + add(a[k]*mu[k],k=1..n))=MRHS;
  #if Debug then  print(newLeq,newMeq) fi;
  Leq:= Leq union {newLeq};
  Meq:= Meq union {newMeq}
od;
Lsoln:= solve(Leq,Lvar);
#if Debug then print(Lsoln) fi;
Msoln:= solve(Meq,Mvar);
L:= lambda[0] + add(lambda[k]*a[k],k=1..n);
M:= mu[0] + add(mu[k]*a[k],k=1..n);
#if Debug then print(Msoln) fi;
if r<s and s<t and t<u 
    then CL:=2+2*esp1(a,n)-2*add(a[k],k=r..u)+add(a[k],k=r..s-1)
              + add(a[k],k=t+1..u-1);
         CMd:= a[u]
elif r<t and t<s and s<u  
    then CL:= 1+esp1(a,n)-add(a[k],k=r..u)+add(a[k],k=t+1..s-1);
         CMd:= 1+esp1(a,n)
elif r<t and t<u and u<s 
     then CL:= 1+ esp1(a,n) - add(a[k],k=r..s-1) + 2*add(a[k],k=t+1..r-1)
           +add(a[k],k=t+1..u-1) + 2*add(a[k],k=u+1..s-1) ;
         CMd:=a[u]
elif t<r and r<s and s<u    
     then CL:= 1+ esp1(a,n) - add(a[k],k=t..u) + add(a[k],k=r..s-1)
      + 2*add(a[k],k=t+1..r-1) ;
          CMd:=a[t]
elif t<r and r<u and u<s
     then CL:= add(a[k],k=t+1..r-1) + add(a[k],k=u+1..s-1);
          CMd:=1+esp1(a,n)
elif t<u and u<r and r<s
     then CL:= add(a[k],k=r..s-1) + add(a[k],k=t+1..u-1)+2*add(a[k],k=u+1..r-1);
          CMd:=a[u]
fi; 
if Verbose then
  print(`The coefficient of`, x[r]*x[s]/x[t]/x[u],` in the n=`,n);
  print(`case of the q-Dyson product`);
  #print(qDysonProd(x,a));
  print(` is conjectured to be`);
  print(
  ``(q,q)[esp1(a,n)]/product(``(q,q)[a[i]],i=1..n) * (1-q^(a[t])) *
   (1-q^(a[u]))
   /(1-q^(1+esp1(a,n)-a[t]))/(1-q^(1+esp1(a,n)-a[u]))/(1-q^(1+esp1(a,n)-a[t]-a[u]))
   *( q^(subs(Lsoln,L)) * (1-q^(1+esp1(a,n) ))
   + q^(subs(Msoln,M)) * (1-q^(1+esp1(a,n)-a[t]-a[u]))    )
   )
else print([r,s,t,u],[subs(Lsoln,L), subs(Msoln,M)-subs(Lsoln,L)],`...`,
            [CL-subs(Lsoln,L), CMd-subs(Msoln,M)+subs(Lsoln,L)]
                       )
fi;
if Debug or Verbose then print(`==========================`) fi;
end: