######################################################################
## BCMV.txt Save this file as  BCMV.txt  to use it,                  #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read  BCMV.txt <Enter>                            #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  DoronZeil at gmail dot com                                       # 
######################################################################


print(`First Written: Feb. 15, 2024: tested for Maple 2020 `):
print(`Version  : Feb. 28, 2024; Adding fascinating code by Mark van Hoeij  `):
print():
print(`This is BCMV.txt, A Maple package`):
print(`accompanying Shalosh B. Ekhad and Doron Zeilberger's  article: `):
print(` Efficient Evaluations of Weighted Sums over the Boolean Lattice inspired by  conjectures of  Berti, Corsi, Maspero, and  Ventura `):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/mamarim/mamarimhtml/bcmv.html `):
print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/BCMV.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():
print(`For general help, and a list of the MAIN functions,`):

print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezra1();`):
print():

print(`For the general procedure, type, ezraG();`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezra1();`):
print():
 



ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are: IV1, IV, IVn, gpJ, npJ `):
print(``):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` BCMV.txt: A Maple package for  investigating a conjectured combinatoiral identity that came up in deep water`):

print(`The MAIN procedures are: A, Ax, AxR, AxC, C, Cx, CxR,  CxH, Bx, BxC, BxR, Dx, DxH, DxR `):
print(` MvHK, K,  ProveConjAx, ProveConjBx, VerifyConjC `):

print(``):

elif nargs=1 and args[1]=A then
print(`A(p): The left side of eq. (3) in the BMV question (as stated). Try: `):
print(`A(10);`):

elif nargs=1 and args[1]=Ax then
print(`Ax(p,x)The left side of eq. (3) in the BMV question where p is replaced by the symbol x. Try:`):
print(`Ax(10,x);`):


elif nargs=1 and args[1]=AxC then
print(`AxC(p,x)The Conjectured Am3x(p,x); Try:`):
print(`AxC(p,x);`):

elif nargs=1 and args[1]=AxR then
print(`AxR(p,x): Ax(p,x) done recursively, rather than directly, via BxR(p,x) (q.v.) . Try:`):
print(`AxR(10,x);`):


elif nargs=1 and args[1]=C then
print(`C(p): The left side of eq. (4) in Alberto Maspero's question. Try:`):
print(`C(10);`):

elif nargs=1 and args[1]=Cx then
print(`Cx(p,x): The left side of eq. (4) in Alberto Maspero's question with p replaced by a symbol x, done directly. Try`):
print(`[seq(Cx(p,x),p=2..10)];`):

elif nargs=1 and args[1]=CxR then
print(`CxR(p,x): The left side of eq. (4) in Alberto Maspero's question with p replaced by a symbol x, done recusively. Same as Cx(p,x) but much faster. `):
print(`[seq(CxR(p,x),p=2..10)];`):

elif nargs=1 and args[1]=CxH then
print(`CxH(p,x): same as CxR(p,x), but faster, using the second-order recurrence `):
print(`[seq(CxH(p,x),p=2..10)];`):

elif nargs=1 and args[1]=Bx then
print(`Bx(p,x): The weight-enumerator (according to npJ(x,L)) of IV(p), done directly. Try:`):
print(`Bx(10,x);`):


elif nargs=1 and args[1]=Dx then
print(`Dx(p,x): The weight-enumerator (according to npJ(x,L)*gpJ(x,L)) of IV(p), done directly. Try:`):
print(`[seq(Dx(i,x),i=1..20)];`):


elif nargs=1 and args[1]=DxH then
print(`DxH(p,x): The weight-enumerator (according to npJ(x,L)*gpJ(x,L)) of IV(p), via the second-order linear recurrence. Same output as Dx(p,x) and DxR(p,x), and  even faster than the latter. Try:`):
print(`[seq(DxH(i,x),i=1..20)];`):

elif nargs=1 and args[1]=DxR then
print(`DxR(p,x): The weight-enumerator (according to npJ(x,L)*gpJ(x,L)) of IV(p), done recursively. Same output as Dx(p,x) but much faster.  Try:`):
print(`[seq(DxR(i,x),i=1..20)];`):

elif nargs=1 and args[1]=BxC then
print(`BxC(p,x): The Conjectured expression for Bx(p,x) Try:`):
print(`BxC(10,x);`):

elif nargs=1 and args[1]=BxR then
print(`BxR(p,x): Bx(p,x) done recursivley (hence much faster). Try:`):
print(`BxR(10,x);`):

elif nargs=1 and args[1]=IV then
print(`IV(p): The set of all increasing vectors of positive integers that end in p. Try:`):
print(`IV(6); `):

elif nargs=1 and args[1]=IV1 then
print(`IV1(p,q): The set of all increasing vectors of length q that end in p. Try:`):
print(`IV1(6,3);`):

elif nargs=1 and args[1]=IVn then
print(`IVn(p): The set of all increasing vectors that end in <p Try:`):
print(`IVn(6);`):

elif nargs=1 and args[1]=gpJ then
print(`gpJ(p,J):inputs an increasing list of pos. integers J and an integr (or symbol p) outputs g_q^(p)(J) in Alberto Maspero's identity. Try:`):
print(`gpJ(11,[3,6,8]);`):


elif nargs=1 and args[1]=K then
 print(`K(n,x): same as CxH(n+1,x)`):

elif nargs=1 and args[1]=MvHK then
print(` MvHK satisfies the same recurrence and initial conditions as K.`):
print(` The first term in the program contains a summation i = 1 .. n`):
print(` but also contains a factor that vanishes at x = n+1.`):
print(` So at x = n+1 only a summation-free formula remains.`):
print(`Try: `):
print(`normal([seq(MvHK(i,x)-K(i,x),i=1..30)]);`):
print(` {seq(MvHK(n-1,n)-n*(n+1)^2,n=2..30)};`):

elif nargs=1 and args[1]=npJ then
print(`npJ(x,J):inputs an increasing list of pos. integers J and an integer (or symbol p) outputs n_q^(p)(J) in the BMV identity. Try:`):
print(`npJ(11,[3,6,8]);`):

elif nargs=1 and args[1]=ProveConjAx then
print(`ProveConjAx(): proves the conjectured expression, AxC(p,x) for Ax(p,x). Try:`):
print(`ProveConjAx() ;`):

elif nargs=1 and args[1]=ProveConjBx then
print(`ProveConjBx(): proves the conjectured expression, BxC(p,x) for Bx(p,x). Try:`):
print(`ProveConjBx() ;`):


elif nargs=1 and args[1]=VerifyConjC then
print(`VerifyConjC(K): verifies the  conjectured expression, p*(p+1)^2, for  Cx(p,p) for p from 2 to K, using CxH(p,p). Try:`):
print(`VerifyConjC(40) ;`):

else
 print(`There is no such thing as`, args):

fi:


end:




#From GuessPol
#GP1(L,n,d): Inputs a list of numbers L and a variable n and a non-neg. d outputs the polynomial of degree d
#P(n), such that P(i)=L[i], for all i from 1 to nops(L). OR FAIL
#nops(L)>=d+2
GP1:=proc(L,n,d) local a,i,P,eq,var:
if nops(L)<d+2 then
 ERROR(`List too short`):
 fi:
#P=a[0]+a[1]*n+...+a[d]*n^d
P:=add(a[i]*n^i,i=0..d):
var:={seq(a[i],i=0..d)}:
eq:={seq(L[i]=subs(n=i,P),i=1..d+2)}:
var:=solve(eq,var):
 if var=NULL then
   RETURN(FAIL):
 fi:
P:=subs(var,P):
if {seq(L[i]-subs(n=i,P),i=d+3..nops(L))}<>{0} then
  RETURN(FAIL):
 fi:
P:
end:



#GP(L,n): Inputs a list of numbers L and a variable n and outputs the polynomial of degree <=nops(L)-2
#P(n), such that P(i)=L[i], for all i from 1 to nops(L). OR FAIL
#nops(L)>=d+2
GP:=proc(L,n) local d,P:
for d from 0 to nops(L)-2 do
 P:=GP1(L,n,d):
 if P<>FAIL then
    RETURN(P):
 fi:
od:
FAIL:
end:

#End From GuessPol

#npJOld(p,J):inputs an increasing list of pos. integers J and an integr (or symbol p) outputs n_q^(p)(J) in Alberto Maspero's identity. Try:
#npJOld(11,[3,6,8]);
npJOld:=proc(p,J) local q,i:
q:=nops(J):
factor((-12)^q*convert(J,`*`)*J[1]*mul(J[i+1]-J[i],i=1..q-1)*(p-J[q])/mul(p^3-p+J[i]-J[i]^3,i=1..q)):
end:



#npJ(x,J):inputs an increasing list of pos. integers J and an integer (or symbol p) outputs n_q^(p)(J) in the BMV identity. Try:
#npJ(11,[3,6,8]);
npJ:=proc(x,J) local q,J1:
option remember:
q:=nops(J):

if q=1 then
RETURN(-12*J[1]^2*(x-J[1])/(x^3-x+J[1]-J[1]^3)):
fi:
J1:=[op(1..q-1,J)]:

factor(normal(npJ(x,J1)*(-12)*J[q]*(J[q]-J[q-1])*(x-J[q])/(x-J[q-1])/(x^3-x+J[q]-J[q]^3))):

end:



#gpJ(x,J):inputs an increasing list of pos. integers J and an integr (or symbol p) outputs g_q^(p)(J) in Alberto Maspero's identity. Try:
#gpJ(11,[3,6,8]);
gpJ:=proc(x,J) local q:

q:=nops(J):
if q=1 then
RETURN(-28/9*x^2+29/45+274/45*J[1]^2-1/6*(x^3-x)/J[1]+1/5*(x^3-x)*(x+J[1])/(x^2+x*J[1]+J[1]^2-1)-13/9*x*J[1]):
fi:
gpJ(x,[op(1..q-1,J)])+
5/18*(x^3-x)*1/J[q]+
38/15*J[q]^2+
(x^3-x)/5*(x+J[q])/(x^2+J[q]^2+x*J[q]-1)
-13/9*J[q-1]*J[q]
-13/9*(J[q]-J[q-1])*x
+49/45
end:

#gpJOld(p,J):inputs an increasing list of pos. integers J and an integr (or symbol p) outputs g_q^(p)(J) in Alberto Maspero's identity. Try:
#gpJOld(11,[3,6,8]);
gpJOld:=proc(x,J) local q,i:
q:=nops(J):

-28/9*x^2+49/45*q+32/9*J[1]^2-4/9-4/9*(x^3-x)/J[1]+ 5/18*(x^3-x)*add(1/J[i],i=1..q)+
38/15*add(J[i]^2,i=1..q)+ (x^3-x)/5*
add((x+J[i])/(x^2+J[i]^2+x*J[i]-1),i=1..q)-13/9*(add(J[i]*J[i+1],i=1..q-1)+J[q]*x):

end:


#IV1(p,q): The set of all increasing vectors of length q that end in p. Try:
#IV1(6,3);
IV1:=proc(p,q) local gu,gu1,p1:
option remember:

if q=1 then
 RETURN({[p]}):
fi:

gu:= {seq(op(IV1(p1,q-1)),p1=1..p-1)}:

gu:={seq([op(gu1),p],gu1 in gu)}:

gu:
end:


#IV(p): The set of all increasing vectors of positive integers that end in p. Try:
#IV(6);
IV:=proc(p) local q:
{seq(op(IV1(p,q)),q=1..p)}:
end:


#IVn(p): The set of all increasing vectors that end in <p Try:
#IVn(6);
IVn:=proc(p) local i:
{seq(op(IV(i)),i=1..p-1)}:
end:


#A(p): The left side of eq. (3) in the BMV question (as stated). Try: A(10);
A:=proc(p) local gu,gu1:
gu:=IV(p) minus {[p]}:
normal(add(npJ(p,[op(1..nops(gu1)-1,gu1)]),gu1 in gu)):
end:




#AxOld(p,x)The left side of eq. (3) in the BMV question where p is replaced by the symbol x
AxOld:=proc(p,x) local gu,gu1:
gu:=IVn(p):
factor(add(npJ(x,gu1),gu1 in gu)):
end:

#Ax(p,x)The left side of eq. (3) in the BMV question where p is replaced by the symbol x
Ax:=proc(p,x) local gu,gu1:
normal(add(Bx(i,x),i=1..p-1)):
end:


#AxC(p,x)The Conjectured Am3x(p,x);
AxC:=proc(p,x):factor((3*p^2*(p-1)^2-2*(2*p-1)*p*(p-1)*x)/x/(x^2-1)):end:

#C(p): The left side of eq. (4) in Alberto Maspero's question
C:=proc(p) local gu,gu1:
gu:=IV(p) minus {[p]}:
add(npJ(p,[op(1..nops(gu1)-1,gu1)])*gpJ(p,[op(1..nops(gu1)-1,gu1)]),gu1 in gu):
end:


#CxOld(p,x): The left side of eq. (4) in Alberto Maspero's question with p replaced by a symbol x
CxOld:=proc(p,x) local gu,gu1:
gu:=IV(p) minus {[p]}:
factor(add(npJ(x,[op(1..nops(gu1)-1,gu1)])*gpJ(x,[op(1..nops(gu1)-1,gu1)]),gu1 in gu)):
end:







#Bx(p,x): The weight-enumerator (according to npJ(x,L)) of IV(p), done directly. Try:
#Bx(10,x);
Bx:=proc(p,x) local gu,gu1:
gu:=IV(p):
factor(add(npJ(x,gu1),gu1 in gu)):
end:

#Hal(p,x)
Hal:=proc(p,x) local gu,i,lu:
lu:=-12*p*(x-p)/(x^3-x+p-p^3):
gu:=p:
for i from 1 to p-1 do
gu:=normal(gu+Bx(i,x)*(p-i)/(x-i)):
od:
normal(lu*gu-Bx(p,x)):
end:



#BxR(p,x): Bx(p,x) done recursivley (hence much faster). Try:
#BxR(10,x);
BxR:=proc(p,x) local i:
option rememeber:
if p=0 then
 RETURN(0):
fi:

	
normal(-12*p*(x-p)/(x^3-x+p-p^3)*(p+add(BxR(i,x)*(p-i)/(x-i),i=1..p-1))):

end:

AxR:=proc(p,x) local i:
normal(add(BxR(i,x),i=1..p-1)):
end:


#BxC(p,x): The Conjectured expression for Bx(p,x) Try:
#BxC(10,x);
BxC:=proc(p,x):
normal(12*p^2*(p-x)/x/(x+1)/(x-1)):
end:

#ProveConjBx(): proves the conjectured expression, BxC(p,x) for Bx(p,x). Try:
#ProveConjBx() ;
ProveConjBx:=proc() local x,p,i:
evalb(0=normal(-12*p*(x-p)/(x^3-x+p-p^3)*(p+sum(BxC(i,x)*(p-i)/(x-i),i=1..p-1))-BxC(p,x))):
end:


#ProveConjAx(): proves the conjectured expression, AxC(p,x) for Ax(p,x). Try:
#ProveConjAx() ;
ProveConjAx:=proc() local i,p,x:
evalb(normal(sum(BxC(i,x),i=1..p-1) -AxC(p,x))=0):
end:


#VerifyConjCold(K): empiricall checks conjecture equation (4) of Alberto Maspero
VerifyConjCold:=proc(K) local p:
evalb({seq(CxR(p,p)-p*(p+1)^2,p=2..K)}={0}):
end:


#VerifyConjC(K): empiricall checks conjecture equation (4) of Alberto Maspero using the recurrence
VerifyConjC:=proc(K) local p:
evalb({seq(CxH(p,p)-p*(p+1)^2,p=2..K)}={0}):
end:



#From GuessPol
#GP1(L,n,d): Inputs a list of numbers L and a variable n and a non-neg. d outputs the polynomial of degree d
#P(n), such that P(i)=L[i], for all i from 1 to nops(L). OR FAIL
#nops(L)>=d+2
GP1:=proc(L,n,d) local a,i,P,eq,var:
if nops(L)<d+2 then
 ERROR(`List too short`):
 fi:
#P=a[0]+a[1]*n+...+a[d]*n^d
P:=add(a[i]*n^i,i=0..d):
var:={seq(a[i],i=0..d)}:
eq:={seq(L[i]=subs(n=i,P),i=1..d+2)}:
var:=solve(eq,var):
 if var=NULL then
   RETURN(FAIL):
 fi:
P:=subs(var,P):
if {seq(L[i]-subs(n=i,P),i=d+3..nops(L))}<>{0} then
  RETURN(FAIL):
 fi:
P:
end:



#GP(L,n): Inputs a list of numbers L and a variable n and outputs the polynomial of degree <=nops(L)-2
#P(n), such that P(i)=L[i], for all i from 1 to nops(L). OR FAIL
#nops(L)>=d+2
GP:=proc(L,n) local d,P:
for d from 0 to nops(L)-2 do
 P:=GP1(L,n,d):
 if P<>FAIL then
    RETURN(P):
 fi:
od:
FAIL:
end:

#End From GuessPol

#npJOld(p,J):inputs an increasing list of pos. integers J and an integr (or symbol p) outputs n_q^(p)(J) in Alberto Maspero's identity. Try:
#npJOld(11,[3,6,8]);
npJOld:=proc(p,J) local q,i:
q:=nops(J):
factor((-12)^q*convert(J,`*`)*J[1]*mul(J[i+1]-J[i],i=1..q-1)*(p-J[q])/mul(p^3-p+J[i]-J[i]^3,i=1..q)):
end:



#npJ(x,J):inputs an increasing list of pos. integers J and an integer (or symbol p) outputs n_q^(p)(J) in the BMV identity. Try:
#npJ(11,[3,6,8]);
npJ:=proc(x,J) local q,J1:
option remember:
q:=nops(J):

if q=1 then
RETURN(-12*J[1]^2*(x-J[1])/(x^3-x+J[1]-J[1]^3)):
fi:
J1:=[op(1..q-1,J)]:

factor(normal(npJ(x,J1)*(-12)*J[q]*(J[q]-J[q-1])*(x-J[q])/(x-J[q-1])/(x^3-x+J[q]-J[q]^3))):

end:



#gpJ(x,J):inputs an increasing list of pos. integers J and an integr (or symbol p) outputs g_q^(p)(J) in Alberto Maspero's identity. Try:
#gpJ(11,[3,6,8]);
gpJ:=proc(x,J) local q:

q:=nops(J):
if q=1 then
RETURN(-28/9*x^2+29/45+274/45*J[1]^2-1/6*(x^3-x)/J[1]+1/5*(x^3-x)*(x+J[1])/(x^2+x*J[1]+J[1]^2-1)-13/9*x*J[1]):
fi:
gpJ(x,[op(1..q-1,J)])+
5/18*(x^3-x)*1/J[q]+
38/15*J[q]^2+
(x^3-x)/5*(x+J[q])/(x^2+J[q]^2+x*J[q]-1)
-13/9*J[q-1]*J[q]
-13/9*(J[q]-J[q-1])*x
+49/45
end:

#gpJOld(p,J):inputs an increasing list of pos. integers J and an integr (or symbol p) outputs g_q^(p)(J) in Alberto Maspero's identity. Try:
#gpJOld(11,[3,6,8]);
gpJOld:=proc(x,J) local q,i:
q:=nops(J):

-28/9*x^2+49/45*q+32/9*J[1]^2-4/9-4/9*(x^3-x)/J[1]+ 5/18*(x^3-x)*add(1/J[i],i=1..q)+
38/15*add(J[i]^2,i=1..q)+ (x^3-x)/5*
add((x+J[i])/(x^2+J[i]^2+x*J[i]-1),i=1..q)-13/9*(add(J[i]*J[i+1],i=1..q-1)+J[q]*x):

end:


#IV1(p,q): The set of all increasing vectors of length q that end in p. Try:
#IV1(6,3);
IV1:=proc(p,q) local gu,gu1,p1:
option remember:

if q=1 then
 RETURN({[p]}):
fi:

gu:= {seq(op(IV1(p1,q-1)),p1=1..p-1)}:

gu:={seq([op(gu1),p],gu1 in gu)}:

gu:
end:


#IV(p): The set of all increasing vectors of positive integers that end in p. Try:
#IV(6);
IV:=proc(p) local q:
{seq(op(IV1(p,q)),q=1..p)}:
end:


#IVn(p): The set of all increasing vectors that end in <p Try:
#IVn(6);
IVn:=proc(p) local i:
{seq(op(IV(i)),i=1..p-1)}:
end:


#A(p): The left side of eq. (3) in the BMV question (as stated). Try: A(10);
A:=proc(p) local gu,gu1:
gu:=IV(p) minus {[p]}:
normal(add(npJ(p,[op(1..nops(gu1)-1,gu1)]),gu1 in gu)):
end:




#AxOld(p,x)The left side of eq. (3) in the BMV question where p is replaced by the symbol x
AxOld:=proc(p,x) local gu,gu1:
gu:=IVn(p):
factor(add(npJ(x,gu1),gu1 in gu)):
end:

#Ax(p,x)The left side of eq. (3) in the BMV question where p is replaced by the symbol x
Ax:=proc(p,x) local gu,gu1:
normal(add(Bx(i,x),i=1..p-1)):
end:


#AxC(p,x)The Conjectured Am3x(p,x);
AxC:=proc(p,x):factor((3*p^2*(p-1)^2-2*(2*p-1)*p*(p-1)*x)/x/(x^2-1)):end:

#C(p): The left side of eq. (4) in Alberto Maspero's question
C:=proc(p) local gu,gu1:
gu:=IV(p) minus {[p]}:
add(npJ(p,[op(1..nops(gu1)-1,gu1)])*gpJ(p,[op(1..nops(gu1)-1,gu1)]),gu1 in gu):
end:


#CxOld(p,x): The left side of eq. (4) in Alberto Maspero's question with p replaced by a symbol x
CxOld:=proc(p,x) local gu,gu1:
gu:=IV(p) minus {[p]}:
factor(add(npJ(x,[op(1..nops(gu1)-1,gu1)])*gpJ(x,[op(1..nops(gu1)-1,gu1)]),gu1 in gu)):
end:







#Bx(p,x): The weight-enumerator (according to npJ(x,L)) of IV(p), done directly. Try:
#Bx(10,x);
Bx:=proc(p,x) local gu,gu1:
gu:=IV(p):
factor(add(npJ(x,gu1),gu1 in gu)):
end:

#Hal(p,x)
Hal:=proc(p,x) local gu,i,lu:
lu:=-12*p*(x-p)/(x^3-x+p-p^3):
gu:=p:
for i from 1 to p-1 do
gu:=normal(gu+Bx(i,x)*(p-i)/(x-i)):
od:
normal(lu*gu-Bx(p,x)):
end:



#BxR(p,x): Bx(p,x) done recursivley (hence much faster). Try:
#BxR(10,x);
BxR:=proc(p,x) local i:
option rememeber:
if p=0 then
 RETURN(0):
fi:

	
normal(-12*p*(x-p)/(x^3-x+p-p^3)*(p+add(BxR(i,x)*(p-i)/(x-i),i=1..p-1))):

end:

AxR:=proc(p,x) local i:
normal(add(BxR(i,x),i=1..p-1)):
end:


#BxC(p,x): The Conjectured expression for Bx(p,x) Try:
#BxC(10,x);
BxC:=proc(p,x):
normal(12*p^2*(p-x)/x/(x+1)/(x-1)):
end:

#ProveConjBx(): proves the conjectured expression, BxC(p,x) for Bx(p,x). Try:
#ProveConjBx() ;
ProveConjBx:=proc() local x,p,i:
evalb(0=normal(-12*p*(x-p)/(x^3-x+p-p^3)*(p+sum(BxC(i,x)*(p-i)/(x-i),i=1..p-1))-BxC(p,x))):
end:




#Dx(p,x): The weight-enumerator (according to npJ(x,L)*gpJ(x,L)) of IV(p), done directly. Try:
#Dx(10,x);
Dx:=proc(p,x) local gu,gu1:
gu:=IV(p):
factor(add(npJ(x,gu1)*gpJ(x,gu1),gu1 in gu)):
end:



#DxRold(p,x): recursive definition of  Dx(p,f1,f2,g1,g2,X,Y). Try:
#DxRold(10,x);
#f1=-12*X^2*(x-X)/(x^3-x+X-X^3),
#f2=(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),
#g1=-28/9*x^2+29/45+274/45*X^2-1/6*(x^3-x)/X+1/5*(x^3-x)*(x+X)/(x^2+x*X+X^2-1)-13/9*x*X,
#g2=5/18*(x^3-x)*1/Y+38/15*Y^2+(x^3-x)/5*(x+Y)/(x^2+Y^2+x*Y-1)-13/9*X*Y-13/9*(Y-X)*x+49/45

DxRold:=proc(p,x) local i:
option remember:
if p=1 then
normal(-12*(x-1)/(x^3-x)*(-28/9*x^2+29/45+274/45-1/6*(x^3-x)+1/5*(x^2-1)-13/9*x))
else
normal(
(-12*p^2*(x-p)/(x^3-x+p-p^3))*(-28/9*x^2+29/45+274/45*p^2-1/6*(x^3-x)/p+1/5*(x^3-x)*(x+p)/(x^2+x*p+p^2-1)-13/9*x*p)+
add(DxRold(i,x)*(-12)*p*(p-i)*(x-p)/(x-i)/(x^3-x+p-p^3),i=1..p-1)+
add(BxR(i,x)*((-12)*p*(p-i)*(x-p)/(x-i)/(x^3-x+p-p^3))*(5/18*(x^3-x)*1/p+38/15*p^2+(x^3-x)/5*(x+p)/(x^2+p^2+x*p-1)-13/9*i*p-13/9*(p-i)*x+49/45),i=1..p-1)):
fi:

end:



#DxR(p,x): recursive definition of  Dx(p,f1,f2,g1,g2,X,Y) and using the explicit expression for Bx(p,x) given by BxC(p,x)
#Try:
#DxR(5,x);
DxR:=proc(p,x) local i,mu,a:
option remember:

mu:=normal(12*a^2*(a-x)/x/(x+1)/(x-1)):

if p=1 then
normal(-12*(x-1)/(x^3-x)*(-28/9*x^2+29/45+274/45-1/6*(x^3-x)+1/5*(x^2-1)-13/9*x))
else
normal(
(-12*p^2*(x-p)/(x^3-x+p-p^3))*(-28/9*x^2+29/45+274/45*p^2-1/6*(x^3-x)/p+1/5*(x^3-x)*(x+p)/(x^2+x*p+p^2-1)-13/9*x*p)+
add(DxR(i,x)*(-12)*p*(p-i)*(x-p)/(x-i)/(x^3-x+p-p^3),i=1..p-1)+
add(subs(a=i,mu)*((-12)*p*(p-i)*(x-p)/(x-i)/(x^3-x+p-p^3))*(5/18*(x^3-x)*1/p+38/15*p^2+(x^3-x)/5*(x+p)/(x^2+p^2+x*p-1)-13/9*i*p-13/9*(p-i)*x+49/45),i=1..p-1)):
fi:

end:

DxH:=proc(n,x) option remember:
if n=1 then
 2/15*(15*x^3+262*x^2+115*x-588)/x/(x+1):
elif n=2 then
4/15/x*(15*x^6+524*x^5+490*x^4-4910*x^3-3649*x^2+1266*x+16056)/(x+1)/(x-1)/(x^2+2*x+3):
else
normal(
2*(100*n^12-26*n^11*x-351*n^9*x^3+78*n^8*x^4+453*n^6*x^6-78*n^5*x^7-199*n^3*x^9+26*n^2*x^10-3*x^12-1200*n^11+286*n^10*x+3159*n^8*x^3-624*n^7*x^4-2718*n^5*x^6+390*n^4*x^7+597*n^2*x^9-52*n*x^10+5900*n^10-897*n^9*x-78*n^8*x^2-11730*n^7*x^3+1122*n^6*x^4+156*n^5*x^5+6642*n^4*x
^6-183*n^3*x^7-78*n^2*x^8-380*n*x^9+12*x^10-15000*n^9-507*n^8*x+624*n^7*x^2+23142*n^6*x^3+2004*n^5*x^4-780*n^4*x^5-8448*n^3*x^6-1011*n^2*x^7+156*n*x^8-18*x^9+19500*n^8+9312*n^7*x-1575*n^6*x^2-26037*n^5*x^3-10086*n^4*x^4+963*n^3*x^5+5655*n^2*x^6+828*n*x^7-18*x^8-7200*n^7-\
23688*n^6*x+714*n^5*x^2+16701*n^4*x^3+15336*n^3*x^4+231*n^2*x^5-1662*n*x^6+54*x^7-13900*n^6+29027*n^5*x+3444*n^4*x^2-5741*n^3*x^3-10868*n^2*x^4-516*n*x^5+12*x^6+21000*n^5-18703*n^4*x-6888*n^3*x^2+771*n^2*x^3+3064*n*x^4-54*x^5-11600*n^4+5784*n^3*x+5265*n^2*x^2+68*n*x^3-3*x^
4+2400*n^3-588*n^2*x-1506*n*x^2+18*x^3)*n/(n-1)/(n-1-x)/(n^2+n*x+x^2-1)/(100*n^9-26*n^8*x-251*n^6*x^3+52*n^5*x^4+202*n^3*x^6-26*n^2*x^7+3*x^9-1350*n^8+312*n^7*x+2259*n^5*x^3-390*n^4*x^4-909*n^2*x^6+78*n*x^7+7800*n^7-1309*n^6*x-52*n^5*x^2-8231*n^4*x^3+740*n^3*x^4+52*n^2*x^5
+1313*n*x^6-61*x^7-25200*n^6+1953*n^5*x+390*n^4*x^2+15501*n^3*x^3+180*n^2*x^4-156*n*x^5-606*x^6+49800*n^5+1601*n^4*x-942*n^3*x^2-15916*n^2*x^3-1482*n*x^4+113*x^5-61650*n^4-9417*n^3*x+729*n^2*x^2+8490*n*x^3+900*x^4+46700*n^3+12874*n^2*x+169*n*x^2-1855*x^3-19800*n^2-7788*n*x
-294*x^2+3600*n+1800*x)*DxH(n-1,x)
-(n^2+n*x+x^2-4*n-2*x+3)*(100*n^9-26*n^8*x-251*n^6*x^3+52*n^5*x^4+202*n^3*x^6-26*n^2*x^7+3*x^9-450*n^8+104*n^7*x+753*n^5*x^3-130*n^4*x^4-303*n^2*x^6+26*n*x^7+600*n^7+147*n^6*x-52*n^5*x^2-701*n^4*x^3-300*n^3*x^4+52*n^2*x^5+101*n*x^6-9*x^7-805*n^5*x+
130*n^4*x^2+147*n^3*x^3+580*n^2*x^4-52*n*x^5-600*n^5+831*n^4*x+98*n^3*x^2+26*n^2*x^3-202*n*x^4+9*x^5+450*n^4-199*n^3*x-277*n^2*x^2+26*n*x^3-100*n^3-52*n^2*x+101*n*x^2-3*x^3)*n/(n-2)/(100*n^9-26*n^8*x-251*n^6*x^3+52*n^5*x^4+202*n^3*x^6-26*n^2*x^7+3*x^9-1350*n^8+312*n^7*x+
2259*n^5*x^3-390*n^4*x^4-909*n^2*x^6+78*n*x^7+7800*n^7-1309*n^6*x-52*n^5*x^2-8231*n^4*x^3+740*n^3*x^4+52*n^2*x^5+1313*n*x^6-61*x^7-25200*n^6+1953*n^5*x+390*n^4*x^2+15501*n^3*x^3+180*n^2*x^4-156*n*x^5-606*x^6+49800*n^5+1601*n^4*x-942*n^3*x^2-15916*n^2*x^3-1482*n*x^4+113*x^5
-61650*n^4-9417*n^3*x+729*n^2*x^2+8490*n*x^3+900*x^4+46700*n^3+12874*n^2*x+169*n*x^2-1855*x^3-19800*n^2-7788*n*x-294*x^2+3600*n+1800*x)/(n^2+n*x+x^2-1)*DxH(n-2,x)):
fi:
end:

#Cx(p,x)The left side of eq. (3) in the BMV question where p is replaced by the symbol x
Cx:=proc(p,x) local i:
normal(add(Dx(i,x),i=1..p-1)):
end:


#CxRold(p,x)The left side of eq. (3) in the BMV question where p is replaced by the symbol x, done recusively
CxRold:=proc(p,x) local i:
normal(add(DxRold(i,x),i=1..p-1)):
end:


#CxR(p,x)The left side of eq. (3) in the BMV question where p is replaced by the symbol x, done recusively
CxR:=proc(p,x) local i:
normal(add(DxR(i,x),i=1..p-1)):
end:



#K(n,x): same as CxH(n+1,x)
K:=proc(n,x) option remember:
if n=1 then
 2/15*(15*x^3+262*x^2+115*x-588)/(x+1)/x:
elif n=2 then
 2/15*(45*x^6+1325*x^5+1372*x^4-10076*x^3-8557*x^2+1599*x+33876)/(x^2+2*x+3)/x/(x+1)/(x-1):
else
normal(
(500*n^13-456*n^12*x-196*n^11*x^2-1706*n^10*x^3+1726*n^9*x^4+946*n^8
*x^5+1402*n^7*x^6-2084*n^6*x^7-1304*n^5*x^8-226*n^4*x^9+814*n^3*x^10+554*n^2*x^
11+30*n*x^12-3250*n^12+2736*n^11*x+1078*n^10*x^2+8530*n^9*x^3-7767*n^8*x^4-3784
*n^7*x^5-4907*n^6*x^6+6252*n^5*x^7+3260*n^4*x^8+452*n^3*x^9-1221*n^2*x^10-554*n
*x^11-15*x^12+6000*n^11-3996*n^10*x-4480*n^9*x^2-12635*n^8*x^3+9184*n^7*x^4+
10758*n^6*x^5+7519*n^5*x^6-6275*n^4*x^7-6810*n^3*x^8-4124*n^2*x^9+547*n*x^10+
262*x^11+2750*n^10-5100*n^9*x+12075*n^8*x^2-640*n^7*x^3+4102*n^6*x^4-19030*n^5*
x^5-6530*n^4*x^6+2130*n^3*x^7+6955*n^2*x^8+3898*n*x^9-70*x^10-19000*n^9+15517*n
^8*x-13586*n^7*x^2+10746*n^6*x^3-14138*n^5*x^4+16200*n^4*x^5+12799*n^3*x^6+7957
*n^2*x^7-2701*n*x^8-1374*x^9+11250*n^8-1372*n^7*x-1253*n^6*x^2+5828*n^5*x^3+
6967*n^4*x^4-5098*n^3*x^5-15122*n^2*x^6-7980*n*x^7+300*x^8+16000*n^7-17938*n^6*
x+13177*n^5*x^2-16440*n^4*x^3-9334*n^3*x^4-6474*n^2*x^5+5499*n*x^6+2550*x^7-\
17750*n^6+7100*n^5*x-7025*n^4*x^2+5360*n^3*x^3+14263*n^2*x^4+6482*n*x^5-330*x^6
-1500*n^5+7589*n^4*x+3227*n^3*x^2+2803*n^2*x^3-5233*n*x^4-2026*x^5+7000*n^4-\
3364*n^3*x-4875*n^2*x^2-1846*n*x^3+115*x^4-2000*n^3-716*n^2*x+1858*n*x^2+588*x^
3)/(n-1)/(n-1-x)/(n^2+n*x+x^2-1)/(250*n^9-228*n^8*x-98*n^7*x^2-603*n^6*x^3+635*
n^5*x^4+375*n^4*x^5+98*n^3*x^6-407*n^2*x^7-277*n*x^8-15*x^9-2250*n^8+1824*n^7*x
+686*n^6*x^2+3618*n^5*x^3-3175*n^4*x^4-1500*n^3*x^5-294*n^2*x^6+814*n*x^7+277*x
^8+8250*n^7-5521*n^6*x-2693*n^5*x^2-8410*n^4*x^3+5715*n^3*x^4+2885*n^2*x^5+946*
n*x^6-362*x^7-15750*n^6+7590*n^5*x+6605*n^4*x^2+9520*n^3*x^3-4445*n^2*x^4-2770*
n*x^5-750*x^6+16500*n^5-4057*n^4*x-9145*n^3*x^2-5870*n^2*x^3+797*n*x^4+965*x^5-\
9000*n^4-324*n^3*x+6503*n^2*x^2+2348*n*x^3+473*x^4+2000*n^3+716*n^2*x-1858*n*x^
2-588*x^3)*K(n-1,x)
-n*(n-x)*(250*n^9-228*n^8*x-98*n^7*x^2-603*n^6*x^3+635*n^5*x^4+375*n
^4*x^5+98*n^3*x^6-407*n^2*x^7-277*n*x^8-15*x^9-750*n^7+863*n^6*x-635*n^5*x^2+
635*n^4*x^3-635*n^3*x^4+635*n^2*x^5+652*n*x^6+45*x^7+750*n^5-1042*n^4*x+635*n^3
*x^2-635*n^2*x^3-473*n*x^4-45*x^5-250*n^3+407*n^2*x+98*n*x^2+15*x^3)*(n^2+n*x+x
^2-2*n-x)/(n-1)/(n-1-x)/(n^2+n*x+x^2-1)/(250*n^9-228*n^8*x-98*n^7*x^2-603*n^6*x
^3+635*n^5*x^4+375*n^4*x^5+98*n^3*x^6-407*n^2*x^7-277*n*x^8-15*x^9-2250*n^8+
1824*n^7*x+686*n^6*x^2+3618*n^5*x^3-3175*n^4*x^4-1500*n^3*x^5-294*n^2*x^6+814*n
*x^7+277*x^8+8250*n^7-5521*n^6*x-2693*n^5*x^2-8410*n^4*x^3+5715*n^3*x^4+2885*n^
2*x^5+946*n*x^6-362*x^7-15750*n^6+7590*n^5*x+6605*n^4*x^2+9520*n^3*x^3-4445*n^2
*x^4-2770*n*x^5-750*x^6+16500*n^5-4057*n^4*x-9145*n^3*x^2-5870*n^2*x^3+797*n*x^
4+965*x^5-9000*n^4-324*n^3*x+6503*n^2*x^2+2348*n*x^3+473*x^4+2000*n^3+716*n^2*x
-1858*n*x^2-588*x^3)*K(n-2,x)):

fi:
end:


# MvHK satisfies the same recurrence and initial conditions as K.
# The first term in the program contains a summation i = 1 .. n
# but also contains a factor that vanishes at x = n+1.
# So at x = n+1 only a summation-free formula remains.
#
MvHK := proc(n,x) local i;
        4*(x-n-1)*(x-n)*(3*n^2+2*n*x+x^2+3*n+x) * add(i/(i^2+i*x+x^2-1)/(x-i), i=1..n)
        +
        n/45*(n+1)*(250*n^4-456*n^3*x-198*n^2*x^2+494*n*x^3-45*x^4+500*n^3-684*n^2*x-198*n*x^2
        +247*x^3+92*n^2-2*n*x+45*x^2-158*n+113*x)/(x^3-x)
end:

#CxH(p,x)The left side of eq. (3) in the BMV question where p is replaced by the symbol x, done via the second-order recurrence
CxH:=proc(p,x)
K(p-1,x):
end:



###start general



ezraG:=proc()
if args=NULL then

print(` BMV.txt: A Maple package for  investigating a conjectured combinatoiral identity that came up in deep water`):

print(`The MAIN procedures are: AxG, AxGr, CxG, BxG, BxGr, DxG, DxGr, HalG, gpJg, npJg `):
print(`  `):

print(``):



elif nargs=1 and args[1]=AxG then
print(`AxG(p,f1,f2,X,Y): Sum of BxG(i,f1,f2,X,Y) for i from 1 to p-1. Try:`):
print(`AxG(10,-12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),X,Y);`):


elif nargs=1 and args[1]=AxGr then
print(`AxGr(p,f1,f2,X,Y): Sum of BxGr(i,f1,f2,X,Y) for i from 1 to p-1. Try:`):
print(`AxGr(10,-12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),X,Y);`):

elif nargs=1 and args[1]=CxG then
print(`CxG(p,f1,f2,g1,g2,X,Y):  add(DxG(i,f1,f2,g1,g2,X,Y),i=1..p-1):. Try:`):
print(`CxG(10, -12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),-28/9*x^2+29/45+274/45*X^2-1/6*(x^3-x)/X+1/5*(x^3-x)*(x+X)/(x^2+x*X+X^2-1)-13/9*x*X,5/18*(x^3-x)*1/Y+38/15*Y^2+(x^3-x)/5*(x+Y)/(x^2+Y^2+x*Y-1)-13/9*X*Y-13/9*(Y-X)*x+49/45,X,Y);`):

elif nargs=1 and args[1]=BxG then
print(`BxG(p,f1,f2,X,Y): The weight-enumerator (according to npJg(J,f1,f2,X,Y)) of IV(p), done directly. Try:`):
print(`BxG(10,-12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),X,Y);`):

elif nargs=1 and args[1]=BxGr then
print(`BxGr(p,f1,f2,X,Y): BxG(p,x) done recursivley (hence much faster). Try:`):
print(`BxGr(10,-12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),X,Y);`):


elif nargs=1 and args[1]=DxG then
print(`DxG(p,f1,f2,g1,g2,X,Y): The weight-enumerator (according to npJg(J,f1,f2,X,Y)*gpJg(J,g1,g2,L)) of IV(p), done directly. Try:`):
print(`DxG(10, -12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),-28/9*x^2+29/45+274/45*X^2-1/6*(x^3-x)/X+1/5*(x^3-x)*(x+X)/(x^2+x*X+X^2-1)-13/9*x*X,5/18*(x^3-x)*1/Y+38/15*Y^2+(x^3-x)/5*(x+Y)/(x^2+Y^2+x*Y-1)-13/9*X*Y-13/9*(Y-X)*x+49/45,X,Y);`):


elif nargs=1 and args[1]=DxGr then
print(`DxGr(p,f1,f2,g1,g2,X,Y): Recusive definition of  DxG(p,f1,f2,g1,g2,X,Y) (q.v.). Try:`):
print(`DxGr(10, -12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),-28/9*x^2+29/45+274/45*X^2-1/6*(x^3-x)/X+1/5*(x^3-x)*(x+X)/(x^2+x*X+X^2-1)-13/9*x*X,5/18*(x^3-x)*1/Y+38/15*Y^2+(x^3-x)/5*(x+Y)/(x^2+Y^2+x*Y-1)-13/9*X*Y-13/9*(Y-X)*x+49/45,X,Y);`):


elif nargs=1 and args[1]=gpJg then
print(`gpJg(J,g1,g2,X,Y):inputs an increasing list of pos. integers J and a function g1 of X and another function g2 of X and Y, outputs`):
print(`the function of J that equals f1(J[1]) of J has length 1, and if J=[J[1],..., J[q]] and J1 is the truncated`):
print(`list [J[1],...,J[q-1]], then it equals its value at J1 plus f2(J[q-1],J[q]). Try:`):
print(`gpJg([3,6,8],-28/9*x^2+29/45+274/45*X^2-1/6*(x^3-x)/X+1/5*(x^3-x)*(x+X)/(x^2+x*X+X^2-1)-13/9*x*X,5/18*(x^3-x)*1/Y+38/15*Y^2+(x^3-x)/5*(x+Y)/(x^2+Y^2+x*Y-1)-13/9*X*Y-13/9*(Y-X)*x+49/45,X,Y);`):

elif nargs=1 and args[1]=HalG then
print(`HalG(p,f1,f2,g1,g2,X,Y); tests the identity that would hopefully lead to a recursive definition of DxG(p,f1,f2,g1,g2). Try:`):
print(`HalG(10, -12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),-28/9*x^2+29/45+274/45*X^2-1/6*(x^3-x)/X+1/5*(x^3-x)*(x+X)/(x^2+x*X+X^2-1)-13/9*x*X,5/18*(x^3-x)*1/Y+38/15*Y^2+(x^3-x)/5*(x+Y)/(x^2+Y^2+x*Y-1)-13/9*X*Y-13/9*(Y-X)*x+49/45,X,Y);`):

elif nargs=1 and args[1]=npJg then
print(`npJg(J,f1,f2,X,Y):inputs an increasing list of pos. integers J functions f1 of X and f2 of X and Y`):
print(`outputs the resulting weight of J, where it defined recursively as follows.`):
print(`if J=[j1] then it is subs(X=j1,f1) and if J=[j[1],..., j[q]] then`):
print(`it npJg([j[1],..., j[q-1])*subs({X=j[q-1],Y=j[q]},f2). Try:`):
print(`npJg([3,6,8], -12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),X,Y);`):





else
 print(`There is no such thing as`, args):

fi:


end:





#npJg(J,f1,f2,X,Y):inputs an increasing list of pos. integers J functions f1 of X and f2 of X and Y
#outputs the resulting weight of J, where it defined recursively as follows.
#if J=[j1] then it is subs(X=j1,f1) and if J=[j[1],..., j[q]] then
#it npJg([j[1],..., j[q-1])*subs({X=j[q-1],Y=j[q]},f2). Try:
#npJg([3,6,8], -12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),X,Y);
npJg:=proc(J,f1,f2,X,Y) local q,J1:
option remember:
q:=nops(J):

if q=1 then
RETURN(factor(subs(X=J[1],f1))):
fi:
J1:=[op(1..q-1,J)]:

factor(normal(npJg(J1,f1,f2,X,Y)*subs({X=J[q-1],Y=J[q]},f2))):

end:




#BxG(p,f1,f2,X,Y): The weight-enumerator (according to npJg(J,f1,f2,X,Y)) of IV(p), done directly. Try:
#BxG(10,-12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),X,Y);
BxG:=proc(p,f1,f2,X,Y) local gu,gu1:
gu:=IV(p):
factor(add(npJg(gu1,f1,f2,X,Y),gu1 in gu)):
end:





#AxG(p,f1,f2,X,Y): Sum of BxG(i,f1,f2,X,Y) for i from 1 to p-1. Try
#AxG(10,-12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),X,Y);
AxG:=proc(p,f1,f2,X,Y) local i:
normal(add(BxG(i,f1,f2,X,Y),i=1..p-1)):
end:





#BxGr(p,f1,f2,X,Y): BxG(p,x) done recursivley (hence much faster). Try:
#BxGr(10,-12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),X,Y);
BxGr:=proc(p,f1,f2,X,Y) local i:
option rememeber:
if p=0 then
 RETURN(0):
fi:

if p=1 then
 RETURN(subs(X=1,f1)):
fi:
	
normal(subs(X=p,f1)+add(BxGr(i,f1,f2,X,Y)*subs({X=i,Y=p},f2),i=1..p-1)):



end:


#AxGr(p,f1,f2,X,Y): Sum of BxGr(i,f1,f2,X,Y) for i from 1 to p-1. Try
#AxGr(10,-12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),X,Y);
AxGr:=proc(p,f1,f2,X,Y) local i:
normal(add(BxGr(i,f1,f2,X,Y),i=1..p-1)):
end:





#gpJg(J,g1,g2,X,Y):inputs an increasing list of pos. integers J and a function g1 of X and another function g2 of X and Y, outputs
#the function of J that equals f1(J[1]) of J has length 1, and if J=[J[1],..., J[q]] and J1 is the truncated
#list [J[1],...,J[q-1]], then it equals its value at J1 plus f2(J[q-1],J[q]). Try:
#gpJg([3,6,8],-28/9*x^2+29/45+274/45*X^2-1/6*(x^3-x)/X+1/5*(x^3-x)*(x+X)/(x^2+x*X+X^2-1)-13/9*x*X,5/18*(x^3-x)*1/Y+38/15*Y^2+(x^3-x)/5*(x+Y)/(x^2+Y^2+x*Y-1)-13/9*X*Y-13/9*(Y-X)*x+49/45,X,Y);
gpJg:=proc(J,g1,g2,X,Y) local q:

q:=nops(J):
if q=1 then
RETURN(subs(X=J[1],g1)):
fi:
gpJg([op(1..q-1,J)],g1,g2,X,Y)+subs({X=J[q-1],Y=J[q]},g2):
end:




#DxG(p,f1,f2,g1,g2,X,Y): The weight-enumerator (according to npJg(J,f1,f2,X,Y)*gpJg(J,g1,g2,L)) of IV(p), done directly. Try:
#DxG(10, -12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),-28/9*x^2+29/45+274/45*X^2-1/6*(x^3-x)/X+1/5*(x^3-x)*(x+X)/(x^2+x*X+X^2-1)-13/9*x*X,5/18*(x^3-x)*1/Y+38/15*Y^2+(x^3-x)/5*(x+Y)/(x^2+Y^2+x*Y-1)-13/9*X*Y-13/9*(Y-X)*x+49/45,X,Y);
DxG:=proc(p,f1,f2,g1,g2,X,Y) local gu,gu1:
gu:=IV(p):
normal(add(npJg(gu1,f1,f2,X,Y)*gpJg(gu1,g1,g2,X,Y),gu1 in gu)):
end:



#CxG(p,f1,f2,g1,g2,X,Y):  add(DxG(i,f1,f2,g1,g2,X,Y),i=1..p-1):. Try:
#CxG(10, -12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),-28/9*x^2+29/45+274/45*X^2-1/6*(x^3-x)/X+1/5*(x^3-x)*(x+X)/(x^2+x*X+X^2-1)-13/9*x*X,5/18*(x^3-x)*1/Y+38/15*Y^2+(x^3-x)/5*(x+Y)/(x^2+Y^2+x*Y-1)-13/9*X*Y-13/9*(Y-X)*x+49/45,X,Y);
CxG:=proc(p,f1,f2,g1,g2,X,Y) local i:
normal(add(DxG(i,f1,f2,g1,g2,X,Y),i=1..p-1)):
end:


#HalG(p,f1,f2,g1,g2,X,Y); tests the identity that would hopefully lead to a recursive definition of DxG(p,f1,f2,g1,g2). Try:
#HalG(10, -12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),-28/9*x^2+29/45+274/45*X^2-1/6*(x^3-x)/X+1/5*(x^3-x)*(x+X)/(x^2+x*X+X^2-1)-13/9*x*X,5/18*(x^3-x)*1/Y+38/15*Y^2+(x^3-x)/5*(x+Y)/(x^2+Y^2+x*Y-1)-13/9*X*Y-13/9*(Y-X)*x+49/45,X,Y);
HalG:=proc(p,f1,f2,g1,g2,X,Y) local i, gu:
gu:=
subs(X=p,f1)*subs(X=p,g1)+
add(DxG(i,f1,f2,g1,g2,X,Y)*subs({X=i,Y=p},f2),i=1..p-1)+
add(BxG(i,f1,f2,X,Y)*subs({X=i,Y=p},f2)*subs({X=i,Y=p},g2),i=1..p-1)-
DxG(p,f1,f2,g1,g2,X,Y):
normal(gu):
end:

#DxGr(p,f1,f2,g1,g2,X,Y): recursive definition of  DxG(p,f1,f2,g1,g2,X,Y). Try:
#DxGr(10, -12*X^2*(x-X)/(x^3-x+X-X^3),(-12)*Y*(Y-X)*(x-Y)/(x-X)/(x^3-x+Y-Y^3),-28/9*x^2+29/45+274/45*X^2-1/6*(x^3-x)/X+1/5*(x^3-x)*(x+X)/(x^2+x*X+X^2-1)-13/9*x*X,5/18*(x^3-x)*1/Y+38/15*Y^2+(x^3-x)/5*(x+Y)/(x^2+Y^2+x*Y-1)-13/9*X*Y-13/9*(Y-X)*x+49/45,X,Y);
DxGr:=proc(p,f1,f2,g1,g2,X,Y) local i:
option remember:
if p=1 then
subs(X=1,f1)*subs(X=1,g1):
else
normal(
subs(X=p,f1)*subs(X=p,g1)
+add(DxGr(i,f1,f2,g1,g2,X,Y)*subs({X=i,Y=p},f2),i=1..p-1)
+add(BxGr(i,f1,f2,X,Y)*subs({X=i,Y=p},f2)*subs({X=i,Y=p},g2),i=1..p-1)
):
fi:

end: