# Written by Aurora Hiveley & Doron Zeilberger, Rutgers University 
# Last updated: 11/1/26 

Help := proc():

if nargs=0 then 

print(`Primary Functions: At(n,t), Bt(n,i,t), AtF(n,t)`):

print(`Statistical Functions: Moms(f,x,r), SMoms(f,x,r)`):

###########################################################################

elif nargs=1 and args[1]=At then 
print(`At(n,t): weight enumerator for the total number of reservoirs across all permutations of length n`):
print(`input: integer n corresponding to permutation length, variable t`):
print(`output: weight enumerator as a polynomial in t for number of reservoirs in all permutations of length n`):
print(`example: At(5,t);`):

elif nargs=1 and args[1]=AtF then 
print(`AtF(n,t): weight enumerator for the total number of reservoirs across all permutations of length n`):
print(`computed according to conjectured holonomic recurrence`):
print(`input: integer n corresponding to permutation length, variable t`):
print(`output: weight enumerator as a polynomial in t for number of reservoirs in all permutations of length n`):
print(`example: AtF(5,t);`):

elif nargs=1 and args[1]=Bt then 
print(`Bt(n,i,t): weight enumerator for the total number of reservoirs across all permutations of length n`):
print(`where element n is in position i`):
print(`input: integer n, integer i, variable t`):
print(`output: weight enumerator as a polynomial in t`):
print(`example: Bt(5,2,t);`):

elif nargs=1 and args[1]=Moms then 
print(`Moms(f,x,r): moments 1 through r of generating function f using variable x`):
print(`input: polynomial f, variable x, integer r (>= 3)`):
print(`output: list of [m1, m2, ..., mr] where mi is the i-th moment of f`):
print(`example: Moms(WE(4,Rpi,t),t,3);`):

elif nargs=1 and args[1]=SMoms then 
print(`SMoms(f,x,r): *standardized* moments 1 through r of generating function f using variable x`):
print(`input: polynomial f, variable x, integer r (>= 3)`):
print(`output: list of [m1, m2, ..., mr] where mi is the i-th standardized moment of f`):
print(`example: SMoms(WE(4,Rpi,t),t,3);`):

fi:

end:

#########################################################################


WhatsNew := proc():
print(`Last updated: 11/1/26`):
print(`Help(); is written and updated`):
end:


#########################################################################


At:=proc(n,t) local i:
if n<=2 then
    RETURN(n!):
fi:

add(Bt(n,i,t),i=1..n):
end:


Bt:=proc(n,i,t) local j:
option remember:

if n=1 or n=2 then
    RETURN(1):
fi:

if i=1 then
    expand(Bt(n-1,1,t)+t*add((j-2)!* binomial(n-2, j-2)*Bt(n-j+1,1,t) ,j=3..n)):

elif i=n then
    expand(t*add((n-j-1)!* binomial(n-2, n-j-1)*Bt(j,j,t),j=1..n-2)+Bt(n-1,n-1,t)):

else
    expand(t*(
    add((abs(i-j)-1)!* binomial(n-2, abs(i-j)-1)*Bt(n-abs(i-j),min(i,j),t) ,j=1..i-2) +
    add((abs(i-j)-1)!* binomial(n-2, abs(i-j)-1)*Bt(n-abs(i-j),min(i,j),t) ,j=i+2..n)   ) 
    + Bt(n-1,i-1,t)+Bt(n-1,i,t) ):
fi:
end:


AtF:=proc(n,t) option remember:
if n<=2 then
    n!:
else
    expand(n*AtF(n-1,t)+ (2*n*t-2*n-4*t+4)*AtF(n-2,t)):
fi:
end:



Moms:=proc(f,x,r) local mu,L,f1,sig,i:

# first moment
f1:=f/subs(x=1,f):
mu:=subs(x=1,diff(f1,x)):
L:=[mu]:

# second moment
f1:=f1/x^mu:
f1:=x*diff(f1,x):
f1:=x*diff(f1,x):
sig:=sqrt(subs(x=1,f1)):
L:=[mu,sig]:

for i from 3 to r do
    f1:=x*diff(f1,x):
    L:=[op(L), subs(x=1,f1)]:
od:

L:
end:


SMoms:=proc(f,x,r) local L,sig,i:
L:=Moms(f,x,r):
sig:=L[2]:
[op(1..2,L),seq(L[i]/sig^i    , i=3..r)]:
end: