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


print(`Written : June 2022  `):
print():
print(`This is  StirlingDet.txt, A Maple package`):
print(`accompanying Tewodros Amdeberhan and Shalosh B. Ekhad's  article: `):
print(` Computing Determinants involving Stirling Numbers  `):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/StirlingDet.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):

print():
print(`------------------------------------`):
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(`For a list of the OBSOLETE functions type: ezraObsolete();`):
print():
 print(`------------------------------------`):




ezraObsolete:=proc()
if args=NULL then

print(`The Obsolete procedures are`):
print(` B1, Bnab,  Bab, babP, nabM, CLD, SEQabSLOW, PaperOld, PaperP `):

else
ezra(args):

fi:


end:

ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are`):
print(`  B1sc, Cgf1, Cgf, CtoR, GuessRec, CtoR, SEQab,  SEQabSLOW, SLx`):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(`  Stirling1Det.txt: A Maple package for  generating and proving explicit experssions for determinants involving (unsigned) Stirling numbers of the first kind`):

print(`The MAIN procedures are`):
print(` B1sc, B1scN, Paper1  `):




elif nargs=1 and args[1]=B1 then
print(`B1(n,a,b): The determinant of the n by n matrix matrix Stirling1(a+i,b+j) done via Dodgson. Faster then Bnab(n,a,b). Try:`):
print(`B1(5,3,1);`):

elif nargs=1 and args[1]=B1sc then
print(`B1sc(a,b,n): B1(n,a,b) according to the Schur formula, Also works for symbolic n, giving a faster way to compute Bnab(a,b,n,K). Try`):
print(`B1sc(5,2,n);`):


elif nargs=1 and args[1]=B1scN then
print(`B1scN(a,b,n): an alternative way to compute B1sc(a,b,n). Try:`):
print(`B1scN(5,2,n);`):

elif nargs=1 and args[1]=B1scM then
print(`B1scM(a,b,n): the matrix on the top of B1sc(a,b,n) `):
print(`B1scM(5,2,n);`):

elif nargs=1 and args[1]=Bab then
print(`Bab(a,b,n,K): Inputs positive integers a and b, a>=b, and outputs a guessed explicit formula for beta_n(a,b)=det(M_n(a,b)) where M_n(a,b) is the determinant of the matrix BnabM(n,a,b)`):
print(`K is the number of data points used. OBSOLETE BECAUSE OF B1s(a,b,n). Try`):
print(`Bab(5,2,n,50);`):


elif nargs=1 and args[1]=BabP then
print(`BabP(a,b,n,K): Inputs positive integers a and b, a>=b, and outputs a guessed explicit formula, and PROOF, for beta_n(a,b)=det(M_n(a,b)) where M_n(a,b) is the determinant of the matrix BnabM(n,a,b)`):
print(`K is the number of data points used. Try`):
print(`BabP(5,2,n,50);`):

elif nargs=1 and args[1]=Bnab then
print(`Bnab(n,a,b): The determinant of the n by n matrix matrix Stirling1(a+i,b+j) done directly. For checking purposes only. Try:`):
print(`Bnab(5,3,1);`):


elif nargs=1 and args[1]=BnabM then
print(`BnabM(n,a,b): The  n by n matrix matrix Stirling1(a+i,b+j) done directly. For checking purposes only. Try:`):
print(`BnabM(5,3,1);`):

elif nargs=1 and args[1]=CLD then
print(`CLD(a,b,K): verfieds the Dodgson identit for Bab(n). Try:`):
print(`CLD(5,2,50);`):


elif nargs=1 and args[1]=Cgf then
print(`Cgf(f,n,q): The generating function for Sum(f*q^n,n=0..infinity) if f is a sum of terms of the form  c*a^n for some numeric c and a`):
print(`Try: Cgf(2^n+3^n,n,q);`):

elif nargs=1 and args[1]=Cgf1 then
print(`Cgf1(f,n,q): The generating function for Sum(f*q^n,n=0..infinity) if f is a constant or of the form c*a^n for some a`):
print(`Try: Cgf1(2^n,n,q);`):


elif nargs=1 and args[1]=CtoR then
print(`CtoR(S,t): Taken from the Maple package Cfinite.txt . `):
print(` The rational function, in t, whose coefficients`):
print(`are the members of the C-finite sequence S. For example, try:`):
print(`CtoR([[1,1],[1,1]],t);`):

elif nargs=1 and args[1]=GuessRec then
print(`GuessRec(L):  Taken from the Maple package Cfinite.txt `):
print(`inputs a sequence L and tries to guess`):
print(`a recurrence operator with constant cofficients `):
print(`satisfying it. It returns the initial  values and the operator`):
print(` as a list. For example try:`):
print(`GuessRec([1,1,1,1,1,1]);`):

elif nargs=1 and args[1]=Paper1 then
print(`Paper1(N,n,q): a paper with explicit expressions for B[a,b](n) for all N>=a>=b>=1. It also gives the generating functions in the variable q. Try:`):
print(`Paper1(5,n,q):`):


elif nargs=1 and args[1]=PaperP then
print(`PaperP(N,K): a paper with explicit expressions for B[a,b](n) for all K>=a>=b>=1. K is the guessing parameter. WITH PROOFS. Try:`):
print(`PaperP(5,100):`):


elif nargs=1 and args[1]=SEQab then
print(`SEQab(a,b,N): the first N+1 terms of beta_n(a,b) starting at n=0. Done using Dogdson, hence faster. Try:`):
print(`SEQab(3,2,100);`):


elif nargs=1 and args[1]=SEQabSLOW then
print(`SEQabSLOW(a,b,N): the first N+1 terms of beta_n(a,b) starting at n=0. Done directly, hence slower. Try:`):
print(`SEQabSLOW(3,2,100);`):

elif nargs=1 and args[1]=SLx then
print(`SL(L,x): Inputs a partition L, and a list of number x, Outouts the Schur polynomial. Try`):
print(`SL([2,1],[x1,x2]);`):


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

fi:


end:





with(linalg):
with(combinat):



##START FROM Cfinite.txt




#GuessRec1(L,d): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients of order d
#satisfying it. It returns the initial d values and the operator
# as a list. For example try:
#GuessRec1([1,1,1,1,1,1],1);
GuessRec1:=proc(L,d) local eq,var,a,i,n:

if nops(L)<=2*d+2 then
 print(`The list must be of size >=`, 2*d+3 ):
 RETURN(FAIL):
fi:

var:={seq(a[i],i=1..d)}:

eq:={seq(L[n]-add(a[i]*L[n-i],i=1..d),n=d+1..nops(L))}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
else
 RETURN([[op(1..d,L)],[seq(subs(var,a[i]),i=1..d)]]):
fi:

end:




#GuessRec(L): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients 
#satisfying it. It returns the initial  values and the operator
# as a list. For example try:
#GuessRec([1,1,1,1,1,1]);
GuessRec:=proc(L) local gu,d:

for d from 1 to trunc(nops(L)/2)-2 do
 gu:=GuessRec1(L,d):
 if gu<>FAIL then
   RETURN(gu):
 fi:
od:
FAIL:

end:

 
#SeqFromRec(S,N): Inputs S=[INI,ope]
#where INI is the list of initial conditions, a ope a list of
#size L, say, and a recurrence operator ope, codes a list of
#size L, finds the first N0 terms of the sequence satisfying
#the recurrence f(n)=ope[1]*f(n-1)+...ope[L]*f(n-L).
#For example, for the first 20 Fibonacci numbers,
#try: SeqFromRec([[1,1],[1,1]],20);
SeqFromRec:=proc(S,N) local gu,L,n,i,INI,ope:
INI:=S[1]:ope:=S[2]:
if not type(INI,list) or not type(ope,list) then
 print(`The first two arguments must be lists `):
 RETURN(FAIL):
fi:
L:=nops(INI):
if nops(ope)<>L then
 print(`The first two arguments must be lists of the same size`):
RETURN(FAIL):
fi:

if not type(N,integer) then
 print(`The third argument must be an integer`, L):
 RETURN(FAIL):
fi:

if N<L then
 RETURN([op(1..N,INI)]):
fi:

gu:=INI:

for n from 1 to N-L do
 gu:=[op(gu),expand(add(gu[nops(gu)-i+1]*ope[i],i=1..L))]:
od:

gu:
end:


#CtoR(S,t), the rational function, in t, whose coefficients
#are the members of the C-finite sequence S. For example, try:
#CtoR([[1,1],[1,1]],t);

CtoR:=proc(S,t) local D1,i,N1,L1,f,f1,L:
if not (type(S,list) and  nops(S)=2 and type(S[1],list) and type(S[2],list)
        and nops(S[1])=nops(S[2]) and type(t, symbol) ) then
   print(`Bad input`):
   RETURN(FAIL):
fi:

D1:=1-add(S[2][i]*t^i,i=1..nops(S[2])):
N1:=add(S[1][i]*t^(i-1),i=1..nops(S[1])):
L1:=expand(D1*N1):
L1:=add(coeff(L1,t,i)*t^i,i=0..nops(S[1])-1):
f:=L1/D1:
L:=degree(D1,t)+10:
f1:=taylor(f,t=0,L+1):
if expand([seq(coeff(f1,t,i),i=0..L)])<>expand(SeqFromRec(S,L+1)) then
print([seq(coeff(f1,t,i),i=0..L)],SeqFromRec(S,L+1)):
 RETURN(FAIL):
else
 RETURN(f):
fi:
end:

##End FROM Cfinite.txt

#Bnab(n,a,b): The determinant of the n by n matrix matrix Stirling1(a+i,b+j) done directly. For checking purposes only. Try:`):
#Bnab(5,3,1);
Bnab:=proc(n,a,b) local i,j:

det([seq([seq(abs(stirling1(i+a,j+b)),j=1..n)],i=1..n)]):
end:


#The matrix
BnabM:=proc(n,a,b) local i,j:

[seq([seq(abs(stirling1(i+a,j+b)),j=1..n)],i=1..n)]:
end:



#SEQabSLOW(a,b,N): the first N+1 terms of beta_n(a,b) starting at n=0/
SEQabSLOW:=proc(a,b,N) local n:

[1,seq(Bnab(n,a,b),n=1..N)]:

end:



#SEQab(a,b,N): the first N+1 terms of beta_n(a,b) starting at n=0
SEQab:=proc(a,b,N) local n:

[1,seq(B1(n,a,b),n=1..N)]:

end:

B1:=proc(n,a,b) 
option remember:

if n=0 then
 RETURN(1):

elif n=1 then
 RETURN(abs(stirling1(a+1,b+1))):
elif a<b then
 RETURN(0):
else
 (B1(n-1,a,b)*B1(n-1,a+1,b+1)-B1(n-1,a+1,b)*B1(n-1,a,b+1))/B1(n-2,a+1,b+1):
fi:
end:


#Bab(a,b,n,K): Inputs positive integers a and b, a>=b, and outputs a guessed explicit formula for beta_n(a,b)=det(M_n(a,b)) where M_n(a,b) is the determinant of the matrix BnabM(n,a,b)
Bab:=proc(a,b,n,K) local f,L,C,q,R,i,n1:
option remember:
L:=SEQab(a,b,K):


C:=GuessRec(L):
if C=FAIL then
 RETURN(FAIL):
fi:


f:=CtoR(C,q):


R:=sort([solve(subs(q=1/q,denom(f)),q)]):
f:=add(subs(q=1/R[i],normal((1-q*R[i])*f))*R[i]^n,i=1..nops(R)):
if {seq(subs(n=n1,f)-L[n1+1],n1=1..nops(L)-1)}<>{0} then
 RETURN(FAIL):
fi:

f:

end:



#CLD(a,b,K): verfieds the Dodgson identit for Bab(n)
CLD:=proc(a,b,K) local gu,gu1,gu2,gu3,n,mu:

gu:=Bab(a,b,n,K):

gu1:=Bab(a-1,b-1,n,K):

gu2:=Bab(a,b-1,n,K):

gu3:=Bab(a-1,b,n,K):

if gu=FAIL or gu1=FAIL or gu2=FAIL or gu3=FAIL then
 RETURN(FAIL):
fi:

mu:=simplify(expand(gu1*gu-subs(n=n+1,gu1)*subs(n=n-1,gu)-gu2*gu3)):

mu:

end:



#BabP(a,b,n,K): The explicit expression for beta_n(a,b) together with the proof. Try:
#BabP(4,2,n,100);
BabP:=proc(a,b,n,K) local gu,B,i,j:
 gu:=Bab(a,b,n,K):
 if gu=FAIL then
  RETURN(FAIL):
 fi:

if CLD(a,b,K)<>0 then 
 RETURN(FAIL):
fi:

print(` Theorem `, [a,b] , `:`):
print(`Let `, B[a,b](n), `be the determinant of the n by n matrix `, Stirling1(a+i,b+j), `1<=i,j<=n . We have `):
print(``):
print(B[a,b](n)=gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(B[a,b](n)=gu):
print(``):


print(``):
print(` Proof: By the Dodgson condensation formula, for any a>=b>=0, the sequence B[a,b](n) satisfies the first-order linear inhomogeneous recurrence`):
print(``):
print(`B[a-1,b-1](n)B[a,b](n)- B[a-1,b-1](n+1) B[a,b](n-1) = B[a,b-1](n) B[a-1,b](n) `):
print(``):
print(`Specializing a=`,a, `and b= `, b, `it follows that the sequence `, B[a,b](n), ` satisfies the first-order linear inhomogeneous recurrence`):
print(``):
print(B[a-1,b-1](n)*B[a,b](n)- B[a-1,b-1](n+1)* B[a,b](n-1) = B[a,b-1](n)*B[a-1,b](n) ):
print(``):
print(`By the already proved Theorems`, [a-1,b-1],[a,b-1],[a-1,b], ` this is the same as `):
print(``):
print(Bab(a-1,b-1,n,K)*B[a,b](n)- subs(n=n+1,Bab(a-1,b-1,n,K))* B[a,b](n-1) -Bab(a,b-1,n,K)*Bab(a-1,b,n,K)=0 ):
print(``):
print(`This recurrence is also satisfied by the proposed expression`, gu, `( check! )`):
print(``):
print(` Since `, B[a,b](1)=subs(n=1,gu), `check!, the theorem follows by induction on n. QED.`):
print(``):
end:


#PaperOld(N,K): a paper with explicit expressions for B[a,b](n) for all K>=a>=b>=1. K is the guessing parameter. WITHOUT PROOFS. Try:
#PaperOld(5,100):

PaperOld:=proc(N,K) local a,b,gu,t0,i,j,n:

t0:=time():
print(`Explicit expressions for the determinant of the n by n determinant Stirling1(a+i,b+j), 1<=i,j<=n, for all 1<=b<=a<=`, N):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

print(`Let `, B[a,b](n), `be the determinant of the n by n matrix `, Stirling1(a+i,b+j), `1<=i,j<=n . We have `):


for a from 1 to N do
for b from 1 to a do
print(``):
gu:=Bab(a,b,n,K):
print(``):
if gu<>FAIL then
print(``):
print(B[a,b](n)=gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(B[a,b](n)=gu):
print(``):
fi:

od:
od:

print(`--------------------------`):
print(``):
print(`We hope that you enjoyed this paper that took`, time()-t0, `seconds to produce. `):
print(``):
end:




#PaperP(N,K): a paper with explicit expressions for B[a,b](n) for all K>=a>=b>=1. K is the guessing parameter. WITH PROOFS. Try:
#PaperP(5,100):

PaperP:=proc(N,K) local a,b,gu,t0,i,j,n:

t0:=time():
print(`Explicit expressions for the determinant of the n by n determinant Stirling1(a+i,b+j), 1<=i,j<=n, for all 1<=b<=a<=`, N):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

print(`Let `, B[a,b](n), `be the determinant of the n by n matrix `, Stirling1(a+i,b+j), `1<=i,j<=n . We have `):


for a from 1 to N do
for b from 1 to a do

BabP(a,b,n,K):

od:
od:

print(`--------------------------`):
print(``):
print(`We hope that you enjoyed this paper that took`, time()-t0, `seconds to produce. `):
print(``):
end:


#SL(x): Inputs a partition L, and a list of number x, Outouts the Schur polynomial
SLx:=proc(L,x) local i,j,a,L1:
a:=nops(x):

if nops(L)>a then
 RETURN(0):
fi:

if nops(L)<a then
 L1:=[op(L),0$(a-nops(L))]:
else
 L1:=L:
fi:

normal(det([seq([seq(x[i]^(L1[j]+a-j),j=1..a)],i=1..a)])/det([seq([seq(x[i]^(a-j),j=1..a)],i=1..a)]))


end:


#B1sc(a,b,n): Bab(n,a,b) according to the Schur formula
B1sc:=proc(a,b,n)  local i:
expand(simplify(expand(SLx([n$(a-b)],[seq(i,i=1..a)])))):
end:



#B1scM(a,b,n): Inputs a partition L,  Outouts the matrix on the top of expression for B1sc(a,b,n)
B1scM:=proc(a,b,n) local i,j,L1:

L1:=[n$(a-b),0$b]:


[seq([seq(i^(L1[j]+a-j),j=1..a)],i=1..a)]:


end:

#B1csN(a,b,n):
B1scN:=proc(a,b,n) local i:
(-1)^binomial(a,2)*simplify(det(B1scM(a,b,n))/mul(i!,i=1..a-1)):
end:


#Paper1(N,n,q): a paper with explicit expressions for B[a,b](n) for all N>=a>=b>=1. 
#Paper1(5,n,q):

Paper1:=proc(N,n,q) local a,b,gu,t0,i,j,lu,L,G,L1,G1:

t0:=time():
print(`Explicit expressions for the determinant of the n by n determinant Stirling1(a+i,b+j), 1<=i,j<=n, for all 1<=b<=a<=`, N):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

print(`Let `, B[a,b](n), `be the determinant of the n by n matrix `, Stirling1(a+i,b+j), `1<=i,j<=n . We have `):

L:=[]:
G:=[]:

for a from 1 to N do
L1:=[]:
G1:=[]:
for b from 0 to a do
print(``):
gu:=B1sc(a,b,n):
L1:=[op(L1),gu]:
print(``):
print(B[a,b](n)=gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(B[a,b](n)=gu):
print(``):
lu:=Cgf(gu,n,q):
G1:=[op(G1),lu]:
print(`The generating function`, Sum(B[a,b](n)*q^n,n=0..infinity),  ` is `):
print(``):
print(lu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(lu):



od:
L:=[op(L),L1]:
G:=[op(G),G1]:
od:

print(`--------------------------`):
print(`To sum up the list of lists of length`, N, `such that its [a,b] entry is B[a,b](n) is `):
print(``):
lprint(L):
print(``):
print(`Also the  list of lists of length`, N, `such that its [a,b] entry is the genearting function of B[a,b](n) w.r.t. q is `):
print(``):
lprint(G):


print(``):
print(`We hope that you enjoyed this paper that took`, time()-t0, `seconds to produce. `):
print(``):
end:





#Cgf1(f,n,q): The generating function for Sum(f*q^n,n=0..infinity) if f is a constant or of the form c*a^n for some a
#Try: Cgf1(2^n,n,q);
Cgf1:=proc(f,n,q) local c,f1,p:

if type(f,numeric) then
 RETURN(f/(1-q)):
fi:


if type(f,`^`)  and op(2,f)=n then
 RETURN(1/(1-op(1,f)*q)):
fi:

if type(f,`^`)  and type(op(2,f),numeric) then
 p:=op(2,f):
 f1:=op(1,f):
if type(f1,`^`)  and op(2,f1)=n then
 RETURN(1/(1-op(1,f1)^p*q)):
else
 RETURN(FAIL):
fi:
fi:


if type(f,`*`) then
 c:=op(1,f):
 if not type(c,numeric) then
   RETURN(FAIL):
 else
RETURN(c*Cgf1(normal(f/c),n,q)):
 fi:
fi:



FAIL:
end:


#Cgf(f,n,q): The generating function for Sum(f*q^n,n=0..infinity) where f is a sum of terms of the form c*a^n
#Try: Cgf(5*2^n+7*5^n,n,q);
Cgf:=proc(f,n,q) local i,f1,gu1,gu,lu:

if not type(f,`+`) then
 RETURN (Cgf1(f,n,q)):
fi:

gu:=0:

for i from 1 to nops(f) do
 f1:=op(i,f):

 gu1:=Cgf1(f1,n,q):

 if gu1=FAIL then
   RETURN(FAIL):
 else
 gu:=gu+gu1:
 fi:

od:

gu:=normal(gu):


lu:=taylor(gu,q=0,101):


if {seq(coeff(lu,q,i)-subs(n=i,f),i=0..100)}<>{0} then
print({seq(coeff(lu,q,i)-subs(n=i,f),i=0..100)}):
 print(gu , `did not work out`):
 RETURN(FAIL):
fi:

gu:

end: