######################################################################
##NCFPS: Save this file as  NCFPS                                    #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read NCFPS<Enter>                                  #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: March 14, 2012
 
print(`Created:  March 14, 2012`):
print(` This is NCFPS `):
print(`to study Non-Commutative Formal Power Series`):
print(`It accompanies the article `):
print(` The Reciprocal of 1+ab+aabb+aaabbb+.. for NON-COMMUTING a nd b, Catalan Numbers and`):
print(` and non commutative quadratic equations. `):
print(`by A. Berenstein, V. Retakh, C. Reuternauer, and D. Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):
 print(`For a list of the procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Pashet, sMul,`):
 print(` SubsMono, SubsMonoR `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: Mul, MulR, RetakhSeries, Reci, `):
 print(` SubsPoly, SubsPolyR,  SE `):
 print(`VerifyRetakh `):


elif nops([args])=1 and op(1,[args])=Khaber then
print(`Khaber(F,G): adds F and G`):

elif nops([args])=1 and op(1,[args])=Mul then
print(`Mul(F,G): the product of two non-comm. polynomials`):
print(`Try: Mul({[4,[a,b]],[-3,[b,a]]},{[15,[a,b,a]]});`):

elif nops([args])=1 and op(1,[args])=MulR then
print(`MulR(F,G,r): the product of two non-comm. polynomials`):
print(`up to degree r`):
print(`Try: MulR({[4,[a,b]],[-3,[b,a]]},{[15,[a,b,a]]},2);`):

elif nops([args])=1 and op(1,[args])=Pashet then
print(`Pashet(F): simplifies the polynomial`):
print(`Try(Pashet([[4,[a,b]]],[[-5,[a,b]]]);`):

elif nops([args])=1 and op(1,[args])=Reci then
print(`Reci(F,r): (1+F)^(-1) up to degree r, For example, try:`):
print(`Reci([[1,[a]]],3);`):

elif nops([args])=1 and op(1,[args])=RetakhSeries then
print(`RetakhSeries(a,b,r): the terms up to degree r of the Retakh Series`):
print(`for (1+ab+aabb+aaabbb+...)^(-1). For example, try:`):
print(`RetakhSeries(a,b,4);`):

elif nops([args])=1 and op(1,[args])=SE then
print(`SE(POL,D1,r): Given a polynomial POL in the letters D1 and`):
print(`other letters, finds all the terms up to degree r`):
print(`of the (non-comm.) formal power series solution of the Eq.`):
print(`D1=POL . For example, try:`):
print(`SE([[1,[x]],[1,[x,D1]],[1,[a,D1,b]],[1,[a,D1,b,D1]]],D1,10);`):

elif nops([args])=1 and op(1,[args])=sMul then
print(`sMul(c,F): multiplying F by the constant c`):


elif nops([args])=1 and op(1,[args])=SubsPoly then
print(`SubsPoly(POL,D1,P): Given a Polynomial, POL,`):
print(`(given as a list of monomials)`):
print(`in an alphabet that contains the letter D1, finds`):
print(`and another polynomial P in the remaining letters, finds`):
print(`the polynomial in the remaining letters obtained`):
print(`by substituting D1=P in POL, for example try:`):
print(`SubsPoly([[5,[a,b,D1,b,D1]]],D1,[[1,[a,b]],[-1,[b,a]]]);`):

elif nops([args])=1 and op(1,[args])=SubsPolyR then
print(`SubsPolyR(POL,D1,P,r): Given a Polynomial, POL,`):
print(`(given as a list of monomials)`):
print(`in an alphabet that contains the letter D1, finds`):
print(`and another polynomial P in the remaining letters, finds`):
print(`the polynomial in the remaining letters, up to degree r, obtained`):
print(`by substituting D1=P in POL, for example try:`):
print(`SubsPolyR([[5,[a,b,D1,b,D1]]],D1,[[1,[a,b]],[-1,[b,a]]],10);`):

elif nops([args])=1 and op(1,[args])=SubsMono then
print(`SubsMono(M,D1,P): Given a monomial [c,[m1,m2, ..., mr]]`):
print(`in an alphabet that contains the letter D1, finds`):
print(`and a polynomial P in the remaining letters, find`):
print(`the polynomial in the remaining letters obtained`):
print(`by substituting D1=P, for example try:`):
print(`SubsMono([5,[a,b,D1,b,D1]],D1,[[1,[a,b]],[-1,[b,a]]]);`):

elif nops([args])=1 and op(1,[args])=SubsMonoR then
print(`SubsMonoR(M,D1,P,r): Given a monomial [c,[m1,m2, ..., mr]]`):
print(`in an alphabet that contains the letter D1, finds`):
print(`and a polynomial P in the remaining letters, find`):
print(`the polynomial, up to degree r, in the remaining letters obtained`):
print(`by substituting D1=P, for example try:`):
print(`SubsMonoR([5,[a,b,D1,b,D1]],D1,[[1,[a,b]],[-1,[b,a]]],10);`):


elif nops([args])=1 and op(1,[args])=VerifyRetakh then
print(`VerifyRetakh(r): empirically verifies Volodia Retakh's`):
print(`conjecture up to r terms. Try:`):
print(`VerifyRetakh(6);`):

else
print(`There is no ezra for`,args):
fi:
 
end:



#sMul(c,F): multiplying F by the constant c
sMul:=proc(c,F) local i:
Pashet([seq([c*F[i][1],F[i][2]],i=1..nops(F))]):
end:

#Pashet(F): simplifies the polynomial
#Try(Pashet([[4,[a,b]]],[[-5,[a,b]]]);
Pashet:=proc(F) local T,gu,f,gu1,mu,mu1,lu,lu1,m,i:
gu:={seq(f[2], f in F)}:

for gu1 in gu do
T[gu1]:=0:
od:

for f in F do
T[f[2]]:=T[f[2]]+f[1]:
od:

mu:=[]:

for gu1 in gu do
 if T[gu1]<>0 then
  mu:=[op(mu),[T[gu1],gu1]] :
 fi:
od:

for mu1 in mu do
T[mu1[2]]:=mu1[1]:
od:

lu:=[seq(mu1[2],mu1 in mu)]:
lu:=sort(lu):

gu:=[seq([T[lu1],lu1],lu1 in lu)]:

m:=max(seq(nops(gu1[2]),gu1 in gu)):

for i from 0 to m do
T[i]:=[]:
od:

for gu1 in gu do
T[nops(gu1[2])]:=[op(T[nops(gu1[2])]),gu1]:
od:

for i from 0 to m do
T[i]:=sort(T[i]):
od:

[seq(op(T[i]),i=0..m)]:
end:

#Mul(F,G): the product of two non-comm. polynomials
#Try: Mul({[4,[a,b]],[-3,[b,a]]},{[15,[a,b,a]]});
Mul:=proc(F,G) local f,g:

Pashet([seq(seq([f[1]*g[1],[op(f[2]),op(g[2])]],f in F), g in G)]):

end:


#MulR(F,G,r): the product of two non-comm. polynomials
#up to degree r
#Try: MulR({[4,[a,b]],[-3,[b,a]]},{[15,[a,b,a]]},2);
MulR:=proc(F,G,r) local i,j,gu:

gu:=[]:
for i from 1 to nops(F) do
for j from 1 to nops(G) do
if nops(F[i][2])+nops(G[j][2])<=r then
 gu:=[op(gu),[F[i][1]*G[j][1],[op(F[i][2]),op(G[j][2])]]]:
fi:
od:
od:

Pashet(gu):

end:

#Khaber(F,G): adds F and G
Khaber:=proc(F,G):
Pashet([op(F),op(G)]):
end:


#Reci(F,r): (1+F)^(-1) up to degree r, For example, try:
#Reci([[1,[a]]],3);
Reci:=proc(F,r) local gu,F1,i,mu:
F1:=[]:

for i from 1 to nops(F) do
if nops(F[i][2])<=r then
  F1:=[op(F1),F[i]]:
 fi:
od:


F1:=Pashet(F1):

F1:=sMul(-1,F1):

mu:=F1:
gu:=Khaber([[1,[]]],mu):



for i from 2 to r do
mu:=MulR(F1,mu,r):
gu:=Khaber(mu,gu):
od:

gu:

end:


#RetakhSeries(a,b,r): the terms up to degree r of the Retakh Series
#for (1+ab+aabb+aaabbb+...)^(-1). For example, try:
#RetakhSeries(a,b,4);
RetakhSeries:=proc(a,b,r) local gu,i:
gu:=[seq([1,[a$i,b$i]],i=1..r)]:
Reci(gu,r):
end:


#SubsMono(M,D1,P): Given a monomial [c,[m1,m2, ..., mr]]
#in an alphabet that contains the letter D1, finds
#and a polynomial P in the remaining letters, find
#the polynomial in the remaining letters obtained
#by substituting D1=P, for example try:
#SubsMono([5,[a,b,D1,b,D1]],D1,[[1,[a,b]],[-1,[b,a]]]);
SubsMono:=proc(M,D1,P) local M1,i,mu:
mu:=[[M[1],[]]]:
M1:=M[2]:
for i from 1 to nops(M1) do
 if M1[i]=D1 then
  mu:=Mul(mu,P):
 else
  mu:=Mul(mu,[[1,[M1[i]]]]):
 fi:
od:
mu:
end:




#SubsMonoR(M,D1,P,r): Given a monomial [c,[m1,m2, ..., mr]]
#in an alphabet that contains the letter D1, finds
#and a polynomial P in the remaining letters, find
#the polynomial, up to degree r, in the remaining letters obtained
#by substituting D1=P, for example try:
#SubsMonoR([5,[a,b,D1,b,D1]],D1,[[1,[a,b]],[-1,[b,a]]],3);
SubsMonoR:=proc(M,D1,P,r) local M1,i,mu:
mu:=[[M[1],[]]]:
M1:=M[2]:
for i from 1 to nops(M1) do
 if M1[i]=D1 then
  mu:=MulR(mu,P,r):
 else
  mu:=MulR(mu,[[1,[M1[i]]]],r):
 fi:
od:
mu:
end:





#SubsPoly(POL,D1,P): Given a Polynomial, POL,
#(given as a list of monomials)
#in an alphabet that contains the letter D1, finds
#and another polynomial P in the remaining letters, find
#the polynomial in the remaining letters obtained
#by substituting D1=P in POL, for example try:
#SubsPoly([[5,[a,b,D1,b,D1]]],D1,[[1,[a,b]],[-1,[b,a]]]);
SubsPoly:=proc(POL,D1,P) local i,gu:
gu:=[]:

for i from 1 to nops(POL) do
gu:=Khaber(gu,SubsMono(POL[i],D1,P)):
od:
gu:
end:


#SubsPolyR(POL,D1,P,r): Given a Polynomial, POL,
#(given as a list of monomials)
#in an alphabet that contains the letter D1, finds
#and another polynomial P in the remaining letters, find
#the polynomial (up to degree r) in the remaining letters obtained
#by substituting D1=P in POL, for example try:
#SubsPolyR([[5,[a,b,D1,b,D1]]],D1,[[1,[a,b]],[-1,[b,a]]],10);
SubsPolyR:=proc(POL,D1,P,r) local i,gu:
gu:=[]:

for i from 1 to nops(POL) do
gu:=Khaber(gu,SubsMonoR(POL[i],D1,P,r)):
od:
gu:
end:


#SE(POL,D1,r): Given a polynomial POL in the letters D1 and
#other letters, finds all the terms up to degree r
#of the (non-comm.) formal power series solution of the Eq.
#D1=POL . For example, try:
#SE([[1,[x]],[1,[x,D1]],[1,[a,D1,b]],[1,[a,D1,b,D1]]],D1,10);
SE:=proc(POL,D1,r) local mu,i:

mu:=[]:
for i from 1 to nops(POL) do
 if not member(D1,{op(POL[i][2])}) and nops(POL[i][2])<=r then
  mu:=[op(mu),POL[i]]:
 fi:
od:

for i from 1 to r do
mu:=SubsPolyR(POL,D1,mu,r):
od:

mu:

end:


#VerifyRetakh(r): empirically verifies Volodia Retakh's
#conjecture up to r terms. Try:
#VerifyRetakh(6);
VerifyRetakh:=proc(r) local gu,a,b,x,lu,i,D1:

gu:=SE([[1,[x]],[1,[x,D1]],[1,[a,D1,b]],[1,[a,D1,b,D1]]],D1,r);
gu:=[[1,[]],op(gu)]:

gu:=MulR([[1,[a]]],gu,r):
gu:=MulR(gu,[[1,[b]]],r):

gu:=SubsPolyR(gu,x,[[1,[a,b]], [-1,[b,a]]],r):

lu:=[seq([1,[a$i,b$i]],i=1..r)]:

lu:=Reci(lu,r):

if lu[1]<>[1,[]] then
 print(`Something is wrong!`):
fi:

lu:=[op(2..nops(lu),lu)]:

lu:=sMul(-1,lu):

evalb(lu=gu):

end: