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

with(combinat):

print(`First Written: Oct. 2023: tested for Maple 2020 `):
print(`Version  : Oct. 2023  `):
print():
print(`This is CalabiWilf.txt, A Maple package`):
print(`accompanying Shalsoh B. Ekhad and Doron Zeilberger's  article: `):
print(` Implementing and Experimenting with the Calabi-Wilf algorithm for radnom selection of a subspace over a finite field `):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/CalabiWilf.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():
 




ezraSt:=proc()
if args=NULL then

print(`The Story procedures are: `):
print(` SimuStory1, SimuStory2, SimuStory3, SimuStory4 `):

else
ezra(args):

fi:


end:

ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are: AAlM, AGF, AnkF, AnkR, BnkrR, Di,  EAlM, EGF,  ExMinWt, FI,  GF1, qbin, MinWt, MinWtD, PeIn, PeInD, Span1, Sqnk, Template, Vecs, WtEnu, WtEnuD `):
print(``):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` CalabiWilf.txt: A Maple package for  implememnting and experimenting with the Calabi-Wilf algorithm for random generation of a subspace over a finite field `):

print(`The MAIN procedures are:`):
print(` AAvM, Av1, EAvM,, MinWtStat, MinWtStatD,  NuSM, Rqnk, SSqnk`):




elif nargs=1 and args[1]=AAlM then
print(`AAlM(n,k,q,B,N,K): inputs non-neg. integer q, pos. integres n and k, with 1<=k<=n, a matrix B, and a variable x, gives the`):
print(`APPROXIMATE average, variance, followed by the third up to the N-th (central and SCALED) moment for the r.v.  number of submatrices (consecutive minor) that are identical to B.`):
print(`taken over the sample space of all row-echelon k by n matrices over GF(q).`):
print(`It uses the Calabi-Wilf algorithm for the simulation. Try:`):
print(`AAlM(4,2,2,[[1]],4,1000);`):

elif nargs=1 and args[1]=AAvM then
print(`AAvM(n,k,q,B,N,K): inputs non-neg. integer q, pos. integres n and k, with 1<=k<=n, a matrix B, and a variable x, gives the`):
print(`APPROXIMATE average, variance, followed by the third up to the N-th (central) moment for the r.v.  number of submatrices (consecutive minor) that are identical to B.`):
print(`taken over the sample space of all row-echelon k by n matrices over GF(q).`):
print(`It uses the Calabi-Wilf algorithm for the simulation. Try:`):
print(`AAvM(4,2,2,[[1]],4,1000);`):



elif nargs=1 and args[1]=AGF then
print(`AGF(n,k,q,B,x,K): inputs non-neg. integer q, pos. integres n and k, with 1<=k<=n, a matrix B, and a variable x, gives the`):
print(`APPROXIMATE prob. generating function for the r.v. number of submatrices (consecutive minor) that are identical to x.`):
print(`sampling the Calabi-Wilf algorithm K times. Try:`):
print(`AGF(5,3,7,[[1]],x,1000);`):

elif nargs=1 and args[1]=AnkF then
print(`AnkF(n,k,q): The number of k-dim subpspaces of GF(q)^n using the formula. Try:`):
print(`AnkF(6,3,q);`):

elif nargs=1 and args[1]=AnkR then
print(`AnkR(n,k,q): The number of k-dim subpspaces of GF(q)^n using the recurrence. Try:`):
print(`AnkR(6,3,q);`):


elif nargs=1 and args[1]=Av1 then
print(`Av1(n,k,q): The exact value of the average number of 1s in  k by n row echelon matrix over GF(q). Try:`):
print(`Av1(10,5,3);`):


elif nargs=1 and args[1]=BnkR then
print(`BnkR(n,k,q): The number of k by n row-echelon matrices with one 1 marked. Try:`):
print(`BnkR(6,3,2);`):


elif nargs=1 and args[1]=Di then
print(`Di(a,b)  outputs true with prob. a/b. Try: Di(1,3);`):

elif nargs=1 and args[1]=EAlM then
print(`EAlM(n,k,q,B,N): inputs non-neg. integer q, pos. integres n and k, with 1<=k<=n, a matrix B, and a variable x, gives the`):
print(`EXACT average, variance, and the third up to the N-th (scaled central) moment for the r.v.  number of submatrices (consecutive minor) that are identical to B.`):
print(`taken over the sample space of all row-echelon k by n matrices over GF(q).`):
print(`WARNING: do not make n too large! If the variance is 0 it does not scale them.Try:`):
print(`EAlM(4,2,2,[[1]],4);`):

elif nargs=1 and args[1]=EAvM then
print(`EAvM(n,k,q,B,N): inputs non-neg. integer q, pos. integres n and k, with 1<=k<=n, a matrix B, and a variable x, gives the`):
print(`EXACT average, variance, followed by the third up to the N-th (central) moment for the r.v.  number of submatrices (consecutive minor) that are identical to B.`):
print(`taken over the sample space of all row-echelon k by n matrices over GF(q).`):
print(`WARNING: do not make n too large! Try:`):
print(`EAvM(4,2,2,[[1]],4);`):

elif nargs=1 and args[1]=EGF then
print(`EGF(n,k,q,B,x): inputs non-neg. integer q, pos. integres n and k, with 1<=k<=n, a matrix B, and a variable x, gives the`):
print(`EXACT prob. generating function for the r.v. number of submatrices (consecutive minor) that are identical to x.`):
print(`WARNING: do not make n too large!. Try`):
print(`EGF(4,2,2,[[1]],x);`):

elif nargs=1 and args[1]=ExMinWt then
print(`ExMinWt(n,k,q,K): picks K random k-dim subspaces of GF(q)^n and finds the average minimal weight, followed largest minimal weight, followed by the champion. Try:`):
print(`ExMinWt(100,5,2,1000);`):

elif nargs=1 and args[1]=FI then
print(`FI(T): Given a template T and a pos. integer q fills in the template T where -1 is from 0 to q. Try:`):
print(`FI([[-1,0],[0,-1]],3);`):


elif nargs=1 and args[1]=GF1 then
print(`GF1(n,k,q,x): The weight-enumerator of the   k by n row-echelon matrices according to the number of 1's. `):
print(`same as EGF(n,k,q,[[1]],x) , but much faster, using recurrence. Try:`):
print(`GF1(10,5,2,x);`):

elif nargs=1 and args[1]=MinWt then
print(`MinWt(M,q): The minimum weight of the subspace of GF(q)^n generated by M (where n:=nops(M[1]). Try:`):
print(`M:=Rqnk(2,10,3); MinWt(M,2);`):

elif nargs=1 and args[1]=MinWtD then
print(`MinWtD(M,q): The minimum weight of the DUAL of a subspace of GF(q)^n generated by M (where n:=nops(M[1]). Try:`):
print(`M:=Rqnk(2,10,3); MinWtD(M,2);`):



elif nargs=1 and args[1]=MinWtStat then
print(`MinWtStat(q,n,k,N,K): inputs pos. integer q (prime or prime power), and pos. integers n and k, 1<=k<=n, and a small pos. integer N, and a large positive integer K`):
print(`uses the Calabi-Wilf algorithm to find `):
print(` (i) upper bounds for the smallest weight in a k-dim subspace of GF(q)^n `):
print(` (ii) lower bounds for the highest smallest weight in a k-dim subspace of GF(q)^n`):
print(` (iii) estimated average, variance, and scaled moments up to the N-th.  `):
print(`(iv): A vector space that achives the largest minimal weight. Try:`):
print(` MinWtStat(2,50,2,4,1000); `):


elif nargs=1 and args[1]=MinWtStatD then
print(`MinWtStatD(q,n,k,N,K): inputs pos. integer q (prime or prime power), and pos. integers n and k, 1<=k<=n, and a small pos. integer N, and a large positive integer K`):
print(`uses the Calabi-Wilf algorithm to find `):
print(` (i) upper bounds for the smallest weight in a DUAL of a k-dim subspace of GF(q)^n `):
print(` (ii) lower bounds for the highest smallest weight in a DUAL of k-dim subspace of GF(q)^n`):
print(` (iii) estimated average, variance, and scaled moments up to the N-th of the above duals.`):
print(`(iv): A vector space that achives the largest minimal weight. Try:`):
print(` MinWtStatD(2,50,2,4,1000); `):

elif nargs=1 and args[1]=NuSM then
print(`NuSM(A,B): How many submatrices (minors) that are B are constained (consecutively in the matrix A). Try:`):
print(`NuSM([[1,0,1],[1,1,0],[0,1,1]],[[1,0],[1,1]]);`):


elif nargs=1 and args[1]=qbin then
print(`qbin(q,n,k): the q-binomial coefficient. Try:`):
print(`qbin(q,5,3);`):


elif nargs=1 and args[1]=PeIn then
print(`PeIn(M,q): The perfection index (a perfect code has value 1) for the code. If the minimal weight is even it returns FAIL. Try:`):
print(`M:=Rqnk(2,10,5); PeIn(M,2);`):


elif nargs=1 and args[1]=PeInD then
print(`PeInD(M,q): The perfection index (a perfect code has value 1) for the DUAL code of M. If the minimal weight is even it returns FAIL. Try:`):
print(`M:=Rqnk(2,10,5); PeInD(M,2);`):

elif nargs=1 and args[1]=Rqnk then
print(`Rqnk(q,n,k): a (uniformly at) random n row-echelon matrices over GF(q)^n representing k-dim subspaces. Try:`):
print(`Rqnk(2,10,5);`):


elif nargs=1 and args[1]=SimuStory1 then
print(`SimuStory1(C,q,MIN,MAX,B,N,K): inputs q for GF(q), positive integer MAX, positive integer N and a large positive integer K`):
print(`runs AAlM(C*k,k,q,B,N,K) for k from MIN to MAX. Try:`):
print(`SimuStory1(2,2,20,25,[[1,1],[1,0]],4,1000);`):


elif nargs=1 and args[1]=SimuStory2 then
print(`SimuStory2(C,q,MIN1,MAX1,N,K): inputs q for GF(q), positive integer MAX, positive integer N and a large positive integer K`):
print(`runs AAlM(C*k,k,q,[[1]],N,K) for k from MIN1 to MAX1 and compares it to the corresponding exact values . Try:`):
print(`SimuStory2(2,2,15,20,6,1000);`):


elif nargs=1 and args[1]=SimuStory3 then
print(`SimuStory3(C,q,MIN1,MAX1,K): inputs q for GF(q), positive integer MAX, positive integer N and a large positive integer K`):
print(`estimates the average number of 1s in k by C*k matrices over GF(q) in echelon form for k from MIN1 to MAX1 and compares to the exact value. Try:`):
print(`SimuStory3(2,2,15,20,1000);`):


elif nargs=1 and args[1]=SimuStory4 then
print(`SimuStory4(q,Ns,Ni,Nf,K,N1,N2): inputs`):
print(`(i) q (for GF(q)), (ii) Ns: starting n (ii) Ni: increment of n (iii) Nf: largest n (iv) K: largest k, (v): N1, the number of samples in one simulation, (vi): N2: the largest (scaled) moment you are`):
print(`interested in:`):
print(`Outputs an article about MinWtStat(q,n,k,N2,N1) (q.v.) for n from Ns to Ni with increments of Ni, and k from 1 to K. Try:`):
print(`SimuStory4(2,10,5,40,4,1000,6);`):

elif nargs=1 and args[1]=Span1 then
print(`Span1(M,n,q): inputs a base M for a vector subspace over GF(q)^n, finds all the vectors. Try:`):
print(`M:=Rqnk(2,10,3); Span1(M,10,2);`):

elif nargs=1 and args[1]=Sqnk then
print(`Sqnk(n,k,q): The set of n by k bases in row-echelon forms of the k-dim subspaces of GF(q)^n. Using Recurrence.Try:`):
print(`Sqnk(4,2,2)`):

elif nargs=1 and args[1]=SSqnk then
print(`SSqnk(q,n,k): The set of n by k bases in row-echelon forms of the k-dim subspaces of GF(q)^n. Using Templates.Try:`):
print(`SSqnk(2,4,2)`):

elif nargs=1 and args[1]=Template then
print(`Template(n,A): Given a subset A of {1,...n} forms the template for row-echelon form where the colums index by A form an identity matrix`):
print(`an entry is -1 when it is not committed. Try:`):
print(`Template(4,[1,3]);`):

elif nargs=1 and args[1]=Vecs then
print(`Vecs(q,n): GF(q)^n. Try: Vecs(3,6);`):

elif nargs=1 and args[1]=WtEnu then
print(`WtEnu(M,q,x): The  weight-enumerator according to the variable x, of the subspace of GF(q)^n generated by M (where n:=nops(M[1]). Try:`):
print(`M:=Rqnk(2,10,3); WtEnu(M,2,x);`):

elif nargs=1 and args[1]=WtEnuD then
print(`WtEnuD(M,q,x): The  weight-enumerator according to the variable x, of the DUAL of the subspace of GF(q)^n generated by M (where n:=nops(M[1]). `):
print(`It uses the MacWilliams identity. Try:`):
print(`M:=Rqnk(2,10,3); WtEnuD(M,2,x);`):

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

fi:


end:


##Start from HISTABRUT
#AveAndMoms(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list whose first entry is the average 
#(alias expectation, alias mean)
#followed by the variance, the third moment (about the mean) and
#so on, until the N-th moment (about the mean).
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#AveAndMoms(((1+x)/2)^100,x,4);

AveAndMoms:=proc(f,x,N) local mu,F,memu1,gu,i:
mu:=simplify(subs(x=1,f)):

if mu=0 then
print(f, `is neither a prob. generating function nor can it be made so`):
RETURN(FAIL):
fi:

F:=f/mu:


memu1:=simplify(subs(x=1,diff(F,x))):

gu:=[memu1]:

F:=F/x^memu1:

F:=x*diff(F,x):
for i from 2 to N do
 F:=x*diff(F,x):
 gu:=[op(gu),simplify(subs(x=1,F))]:
od:

gu:

end:


#Alpha(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list, of length N, whose 
#(i) First entry is the average
#(ii): Second entry is the variance
#for i=3...N, the i-th entry is the so-called alpha-coefficients
#that is the i-th moment about the mean divided by the
#variance to the power i/2 (For example, i=3 is the skewness
#and i=4 is the Kurtosis)
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#Alpha(((1+x)/2)^100,x,4);
Alpha:=proc(f,x,N) local gu,i:
 gu:=AveAndMoms(f,x,N):
 if gu=FAIL then
  RETURN(gu):
 fi:

if gu[2]=0 then
#print(`The variance is 0`):
  RETURN(gu):
fi:

[gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..N)]:

end:


##end from HISTABRUT

#Di(a,b)  outputs true with prob. a/b. Try: Di(1,3);
Di:=proc(a,b) evalb(rand(1..b)()<=a):end:

#qbin(q,n,k): the q-binomial coefficient. Try:
#qbin(q,5,3);
qbin:=proc(q,n,k) local i:
normal(mul(1-q^(n-i),i=0..k-1)/mul(1-q^i,i=1..k)):
end:

#Vecs(q,n): GF(q)^n. Try: Vecs(3,6);
Vecs:=proc(q,n) local gu,gu1,i:
option remember:
if n=0 then
 RETURN({[]}):

else
gu:=Vecs(q,n-1):
RETURN({seq(seq([op(gu1),i],gu1 in gu),i=0..q-1)}):
fi:
end:


#Template(n,A): Given a subset A of {1,...n} forms the template for row-echelon form where the colums index by A form an identity matrix
#an entry is -1 when it is not committed. Try:
#Template(4,[1,3]);
Template:=proc(n,A) local i,j,T,k,i1:
k:=nops(A):

for i from 1 to k do
 for j from 1 to n do
  T[i,j]:=-1:
 od:
od:

for i from 1 to k do
 
 for i1 from 1 to k do
  T[i1,A[i]]:=0:
 od:

T[i,A[i]]:=1:

 for j from A[i]+1 to n do
  T[i,j]:=0:
 od:
od:

[seq([seq(T[i,j],j=1..n)],i=1..k)]:

end:

#FI(T): Given a template T and a pos. integer q fills in the template T where -1 is from 0 to q. Try:
#FI([[-1,0],[0,-1]],3);
FI:=proc(T,q) local k,n,i,j,i0,r,T1,gu,j0:
 option remember:

k:=nops(T):
n:=nops(T[1]):

if not member(-1,{seq(seq(T[i,j],j=1..n),i=1..k)}) then
 RETURN({T}):
fi:

for i from 1 to k while not member(-1,{seq(T[i,j],j=1..n)}) do od:

if i=k+1 then
RETURN({T}):
fi:

i0:=i:

for j from 1 to n while T[i0,j]<>-1 do od:
j0:=j:

gu:={}:

for r from 0 to q-1 do
 T1:=[op(1..i0-1,T),[op(1..j0-1,T[i0]),r,op(j0+1..n,T[i0])],op(i0+1..k,T)]:
 gu:=gu union FI(T1,q):
od:
gu:

end:



#SSqnk(q,n,k): All the k by n row-echelon matrices over GF(q)^n representing k-dim subspaces. Try
#SSqnk(2,4,2);
SSqnk:=proc(q,n,k) local i,lu:
lu:=choose([seq(i,i=1..n)],k):
 {seq(op(FI(Template(n,lu[i]),q)),i=1..nops(lu))}:
end:

Sqnk:=proc(q,n,k) local gu1,gu2,gu,lu,i1,gu2a,lu1:
option remember:

if k>n  then
 RETURN({}):
fi:

if k=1 then
 RETURN(SSqnk(q,n,1)): 
fi:

gu1:=Sqnk(q,n-1,k-1):

gu:={seq([[1,0$(n-1)],seq([0,op(gu1[i1])],i1=1..nops(gu1))],gu1 in gu1)}:

gu2:=Sqnk(q,n-1,k):
lu:=Vecs(q,k):

for gu2a in gu2 do
for lu1 in lu do
 gu:=gu union {[seq([lu1[i1],op(gu2a[i1])],i1=1..k)]}:
od:
od:
gu:

end:



#Rqnk(q,n,k): a (uniformly at) random n row-echelon matrices over GF(q)^n representing k-dim subspaces. Try
#Rqnk(2,10,5);
Rqnk:=proc(q,n,k) local ra1,gu1,gu2,gu,lu,i1,i:

if k>n  then
 RETURN({}):
fi:

ra1:=rand(0..q-1):

if k=1 then
 
if n=1 then
  RETURN([[1]]):
  else
   if Di(q-1,q^n-1) then
     RETURN([[1,0$(n-1)]]):
   else
    gu:=Rqnk(q,n-1,1):
   RETURN([[ra1(),op(gu[1])]]):
 fi:

fi:
fi:

if Di(q^k-1,q^n-1) then
 
gu1:=Rqnk(q,n-1,k-1):

gu:=RETURN([[1,0$(n-1)],seq([0,op(gu1[i1])],i1=1..nops(gu1))]):

else
gu2:=Rqnk(q,n-1,k):
lu:=[seq(ra1(),i=1..k)]:

RETURN([seq([lu[i1],op(gu2[i1])],i1=1..k)]):
fi:

end:

#NuSM(A,B): How many submatrices (minors) that are B are constained (consecutively in the matrix A).
#Try:
#NuSM([[1,0,1],[1,1,0],[0,1,1]],[[1,0],[1,1]]);
NuSM:=proc(A,B) local i,j,co,i1,j1:
co:=0:
for i from 1 to nops(A)-nops(B)+1 do
 for j from 1 to nops(A[1])-nops(B[1])+1 do
   if  [seq([seq(A[i1][j1],j1=j..j+nops(B[1])-1)],i1=i..i+nops(B)-1)]=B then
    co:=co+1:
   fi:
 od:
od:
co:
end:


#EGF(n,k,q,B,x): inputs non-neg. integer q, pos. integres n and k, with 1<=k<=n, a matrix B, and a variable x, gives the
#EXACT prob. generating function for the r.v. number of submatrices (consecutive minor) that are identical to x.
#WARNING: do not make n too large!
EGF:=proc(n,k,q,B,x) local gu,gu1:
gu:=SSqnk(q,n,k):
add(x^NuSM(gu1,B),gu1 in gu)/nops(gu):
end:




#EAvM(n,k,q,B,N): inputs non-neg. integer q, pos. integres n and k, with 1<=k<=n, a matrix B, and a variable x, gives the
#EXACT average, variance, and the third up to the N-th (central) moment for the r.v.  number of submatrices (consecutive minor) that are identical to x.
#WARNING: do not make n too large! Try:
#EAvM(2,4,2,[[1]],4);
EAvM:=proc(n,k,q,B,N) local x:
AveAndMoms(EGF(n,k,q,B,x),x,N):
end:


#EAlM(n,k,q,B,N): inputs non-neg. integer q, pos. integres n and k, with 1<=k<=n, a matrix B, and a variable x, gives the
#EXACT average, variance, and the third up to the N-th (scaled central) moment for the r.v.  number of submatrices (consecutive minor) that are identical to x.
#WARNING: do not make n too large! If the variance is 0 it does not scale them.Try:
#EAlM(4,2,5,[[1]],4);
EAlM:=proc(n,k,q,B,N) local x:
Alpha(EGF(n,k,q,B,x),x,N):
end:





#AGF(n,k,q,B,x,K): inputs non-neg. integer q, pos. integres n and k, with 1<=k<=n, a matrix B, and a variable x, gives the
#APPROXIMATE prob. generating function for the r.v. number of submatrices (consecutive minor) that are identical to x.
#sampling the Calabi-Wilf algorithm K times. Try:
#AGF(2,5,3,[[1]],x,1000);
AGF:=proc(n,k,q,B,x,K) local i,M,gu:
gu:=0:
for i from 1 to K do
M:=Rqnk(q,n,k):
gu:=gu+x^NuSM(M,B):
od:
evalf(gu/K):
end:



#AAvM(n,k,q,B,N,K): inputs non-neg. integer q, pos. integres n and k, with 1<=k<=n, a matrix B, and a variable x, gives the
#APPROXIMATE average, variance, and the third up to the N-th (central) moment for the r.v.  number of submatrices (consecutive minor) that are identical to x.
#AAvM(4,2,q,[[1]],4,1000);
AAvM:=proc(n,k,q,B,N,K) local x:
AveAndMoms(AGF(n,k,q,B,x,K),x,N):
end:


#AAlM(n,k,q,B,N,K): inputs non-neg. integer q, pos. integres n and k, with 1<=k<=n, a matrix B, and a variable x, gives the
#APPROXIMATE average, variance, and the third up to the N-th (scaled central) moment for the r.v.  number of submatrices (consecutive minor) that are identical to B
#using simulation K times
#AAlM(2,4,2,[[1]],4,1000);
AAlM:=proc(n,k,q,B,N,K) local x:
Alpha(AGF(n,k,q,B,x,K),x,N):
end:


#SimuStory1(C,q,MIN1,MAX1,B,N,K): inputs q for GF(q), positive integer MAX, positive integer N and a large positive integer K
#runs AAlM(C*k,k,q,B,N,K) for k from MIN1 to MAX1. Try:
#SimuStory1(2,2,15,20,6,[[1]],1000);

SimuStory1:=proc(C,q,MIN1,MAX1,B,N,K) local k,t0:

t0:=time():
print(`Estimating the average, variance, and central scaled moments up to the`, N, `for the number of occurences of the submatrix`, B, ` in row-echelon matrices of dimension k by `,C*k , `matrices over GF(q) with q=`,q ):
print(`for k from`, MIN1,  ` to `, MAX1):
print(`by simulating`, K, `times, using the amazing Calabi-Wilf algorithm for randomly generating a k-subspace of GF(q)^n`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

print(`Below for the given dimension we repeat the simulation using`, K, `samples , three times to check that it is reliable`):

for k from MIN1 to MAX1 do
print(`For `, k, ` by `, C*k , `row-echelon matrices over GF(q) with q=`, q, `the average, variance, and central scaled moments are estimated as follows, doing it three times `):

lprint(AAlM(C*k,k,q,B,N,K) ):
print(``):
lprint(AAlM(C*k,k,q,B,N,K) ):
print(``):
lprint(AAlM(C*k,k,q,B,N,K) ):
print(``):
od:

print(``):
print(`------------------------------------------------------`):
print(``):
print(`This ends this article that took`, time()-t0, `to generate `):
print(``):
end:

#AnkR(n,k,q): The number of k-dim subpspaces of GF(q)^n using the recurrence. Try:
#AnkR(6,3,q);
AnkR:=proc(n,k,q) option remember:
if k=1 then
 RETURN((q^n-1)/(q-1)):
elif  k>n then
 RETURN(0):
else
RETURN(normal(AnkR(n-1,k-1,q)+q^k*AnkR(n-1,k,q))):
fi:
end:


#AnkF(n,k,q): The number of k-dim subpspaces of GF(q)^n using the formula. Try:
#AnkF(6,3,q);
AnkF:=proc(n,k,q) local i:
option remember:
normal(mul(1-q^(n-i),i=0..k-1)/mul(1-q^i,i=1..k)):
end:


#BnkR(n,k,q): The number of k by n row-echelon matrices with one 1 marked. Try
#BnkR(6,3,q);
BnkR:=proc(n,k,q) option remember:
if k=1 then
RETURN((q^n*n*q+q^n*q^2-q^n*n-2*q^n*q-q^2+2*q)/(q-1)^2/q):
# RETURN(n*q^(n-1)):
elif  k>n then
 RETURN(0):
else
RETURN(normal(AnkR(n-1,k-1,q)+BnkR(n-1,k-1,q)+k*q^(k-1)*AnkR(n-1,k,q))+q^k*BnkR(n-1,k,q) ):
fi:
end:

#Av1(n,k,q): The exact value of the average number of 1s in  k by n row echelon matrix over GF(q). Try:
#Av1(n,k,q);
Av1:=proc(n,k,q)
BnkR(n,k,q)/AnkR(n,k,q):
end:




#GF1(n,k,q,x): The weight-enumerator of the   k by n row-echelon matrices according to the number of 1's. 
#GF1(6,3,q,x); 
GF1:=proc(n,k,q,x) option remember:
if k=1 then
 RETURN(normal(x*((x+(q-1))^n-1)/(x+q-2)*(q-1) /(q^n-1)    )    ):
elif  k>n then
 RETURN(0):
else
RETURN(expand(x*GF1(n-1,k-1,q,x)*(q^k-1)/(q^n-1)+(x+(q-1))^k*(q^(n-k)-1)/(q^n-1)*GF1(n-1,k,q,x))):
fi:
end:


#Av1s(n,k,q): The exact value of the average number of 1s in  k by n row echelon matrix over GF(q). Try:
#Av1s(n,k,q);
Av1s:=proc(n,k,q) local gu,x:
gu:=GF1(n,k,q,x):
subs(x=1,diff(gu,x))/subs(x=1,gu):
end:





#SimuStory2(C,q,MIN1,MAX1,N,K): inputs q for GF(q), positive integer MAX, positive integer N and a large positive integer K
#runs AAlM(C*k,k,q,[[1]],N,K) for k from MIN1 to MAX1 and compares it to the corresponding exact values . Try:
#SimuStory2(2,2,15,20,6,1000);

SimuStory2:=proc(C,q,MIN1,MAX1,N,K) local k,t0,f,x:

t0:=time():
print(` Estimating  the average, variance, and central scaled moments up to the`, N, `for the number of occurences of the number of ONES in row-echelon matrices of dimension k by `,C*k , `matrices over GF(q) with q=`,q ):
print(`for k from`, MIN1,  ` to `, MAX1):
print(`by simulating`, K,  `times using the amazing Calabi-Wilf algorithm for randomly generating a k-subspace of GF(q)^n`):
print(`and comparing it to the exact values `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

print(`Below for the given dimension we repeat the simulation, each of them`, K, `times,  three times to check that it is reliable`):

for k from MIN1 to MAX1 do
print(`For `, k, ` by `, C*k , `row-echelon matrices over GF(q) with q=`, q, `the average, variance, and central scaled moments are estimated as follows, doing it three times `):

lprint(AAlM(C*k,k,q,[[1]],N,K) ):
print(``):
lprint(AAlM(C*k,k,q,[[1]],N,K) ):
print(``):
lprint(AAlM(C*k,k,q,[[1]],N,K) ):
print(``):
f:=GF1(C*k,k,q,x):
print(``):
print(`The exact values are `):
print(``):
print(evalf(Alpha(f,x,N)) ):

od:

print(``):
print(`------------------------------------------------------`):
print(``):
print(`This ends this article that took`, time()-t0, `to generate `):
print(``):
end:


#SimuStory3(C,q,MIN1,MAX1,K): inputs q for GF(q), positive integer MAX, positive integer N and a large positive integer K
#estimates the average number of 1s in k by C*k matrices over GF(q) in echelon form for k from MIN1 to MAX1 and compares to the exact value. Try:
#SimuStory3(2,2,15,20,6,1000);

SimuStory3:=proc(C,q,MIN1,MAX1,K) local k,t0,f,x:

t0:=time():
print(`Comparing the Estimates of the average number  of ONES in row-echelon matrices of dimension k by `,C*k , `matrices over GF(q) with q=`,q ):
print(`for k from`, MIN1,  ` to `, MAX1):
print(`by simulating`, K,  `times, using the amazing Calabi-Wilf algorithm for randomly generating a k-subspace of GF(q)^n`):
print(`and comparing it to the exact values `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

print(`Below for the given dimension we repeat the simulations each of them`, K, ` times , and repeating each simulation three times to check that it is reliable`):

for k from MIN1 to MAX1 do
print(`For `, k, ` by `, C*k , `row-echelon matrices over GF(q) with q=`, q, `the average is  estimated as follows, doing it three times `):

lprint(AAvM(C*k,k,q,[[1]],1,K)[1] ):
print(``):
lprint(AAvM(C*k,k,q,[[1]],1,K)[1] ):
print(``):
lprint(AAvM(C*k,k,q,[[1]],1,K)[1] ):
print(``):
f:=evalf(Av1(C*k,k,q)):
print(``):
print(`The exact value is `, f):

od:

print(``):
print(`------------------------------------------------------`):
print(``):
print(`This ends this article that took`, time()-t0, `to generate `):
print(``):
end:

##start Span
#Span1(M,n,q): inputs a base M for a vector subspace over GF(q)^n, finds all the vectors. Try:
#M:=Rqnk(2,10,3); Span1(M,10,2);
Span1:=proc(M,n,q) local k,gu,v,gu1,i,M1:
k:=nops(M):

if k=0 then
 RETURN({[0$n]}):
fi:
v:=M[k]:
M1:=[op(1..k-1,M)]:
gu:=Span1(M1,n,q):
{seq(seq(gu1+i*v mod q, i=0..q-1),gu1 in gu)}:

end:

#Wt1(v): the number of non-zero entries in the vector v
Wt1:=proc(v) local i,co:
co:=0:
for i from 1 to nops(v) do
 if v[i]<>0 then
   co:=co+1:
 fi:
od:
co:
end:

#MinWt(M,q): The minimum weight of the subspace of GF(q)^n generated by M (where n:=nops(M[1]). Try:
#M:=Rqnk(2,10,3); NinWt(M,2);
MinWt:=proc(M,q) local n,gu,v:
n:=nops(M[1]):
gu:=Span1(M,n,q) minus {[0$n]}:
min(seq(Wt1(v),v in gu)):
end:


#ExMinWt(n,k,q,K): picks K random k-dim subspaces of GF(q)^n and finds the average minimal weight, followed largest minimal weight, followed by the champion. Try:
#ExMinWt(100,5,2,1000);
ExMinWt:=proc(n,k,q,K) local aluf,si,su,i,M,hal:
su:=0:
aluf:=Rqnk(q,n,k):
si:=MinWt(aluf,q):
su:=si:
for i from 1 to K-1 do
 M:=Rqnk(q,n,k):
 hal:=MinWt(M,q):
su:=su+hal:
 if hal>si then
   si:=hal:
  aluf:=M:
 fi:
od:

[evalf(su/K),si,aluf]:
end:


#MinWtStat(q,n,k,N,K): inputs pos. integer q (prime or prime power), and pos. integers n and k, 1<=k<=n, and a small pos. integer N, and a large positive integer K
#uses the Calabi-Wilf algorithm to find 
#(i) upper bounds for the smallest weight in a k-dim subspace of GF(q)^n
#(ii) lower bounds for the highest smallest weight in a k-dim subspace of GF(q)^n
#(iii) estimated average, variance, and scaled moments up to the N-th. Try:
#(iv): A vector space that achives the largest minimal weight
#MinWtStat(2,50,2,4,1000);
MinWtStat:=proc(q,n,k,N,K) local f,x,M,i,lu,f1,si,aluf:
f:=0:
si:=1:
aluf:=0:
for i from 1 to K do
 M:=Rqnk(q,n,k):
 lu:=MinWt(M,q):
if lu>si then
  si:=lu:
 aluf:=M:
fi:
 f:=f+x^lu:
od:

f1:=f/K:

[[ldegree(f,x),degree(f,x),evalf(Alpha(f1,x,N))],aluf]:
end:




#SimuStory4(q,Ns,Ni,Nf,K,N1,N2): inputs
#(i) q (for GF(q)), (ii) Ns: starting n (ii) Ni: increment of n (iii) Nf: largest n (iv) K: largest k, (v): N1, the number of samples in one simulation, (vi): N2: the largest (scaled) moment you are
#interested in:
#Outputs an article about MinWtStat(q,n,k,N2,N1) (q.v.) for n from Ns to Ni with increments of Ni, and k from 1 to K. Try:
#SimuStory4(2,10,5,40,4,1000,6);
SimuStory4:=proc(q,Ns,Ni,Nf,K,N1,N2) local n,k,t0,lu:

t0:=time():
print(`On the Minimal Weight of vectors in k-dimensional subspaces of GF(q)^n for k from 1 to`, K, ` q=`, q, ` and n from `, Ns, ` to `, Nf, `in increments of`, Ni):
print(`By Shalosh B. Ekhad `):
print(``):

print(`This article reports simulating the Calabi-Wilf algorithm`, N1, `times for investigating the minimum weight of vectors in low-dimensional subspaces of GF(q)^n`):
print(`with q=`, q):
print(`We run each simulation twice, and for each report the smallest smallest weight recorded in the sample (offering upper bound for the smallest possible smallest weight`):
print(`followed by the largest smallest weight, offering lower bounds for this important quantity`):
print(`followed by the estimated average, variance, and scaled central moments up to the`,N2):
print(``):

for n from Ns by Ni to Nf do
print(`Investigating n=`, n):
print(``):
for k from 1 to K do
print(``):
print(`For `, k, `dim. subspaces of GF(q)^n  with q=`, q,  ` with n= `, n, `the numbers are:`):
print(``):
lu:=MinWtStat(q,n,k,N2,N1):
lprint(lu[1]):
print(``):
lu:=MinWtStat(q,n,k,N2,N1):
print(`and again:`):
print(``):
lprint(lu[1]):
print(``):
od:
od:

print(`------------------------------------------------------`):
print(``):
print(`This ends this article that took`, time()-t0, `to generate `):
print(``):
end:






#WtEnu(M,q,x): The  weight-enumerator according to the variable x, of the subspace of GF(q)^n generated by M (where n:=nops(M[1]). Try:
#M:=Rqnk(2,10,3); WtEnu(M,2,x);
WtEnu:=proc(M,q,x) local n,gu,v:
n:=nops(M[1]):
gu:=Span1(M,n,q) minus {[0$n]}:
1+add(x^Wt1(v),v in gu):
end:


#WtEnuD(M,q,x): The  weight-enumerator according to the variable x, of the DUAL of the subspace of GF(q)^n generated by M (where n:=nops(M[1]). 
#It uses the MacWilliams identity Try:
#M:=Rqnk(2,10,3); WtEnuD(M,2,x);
WtEnuD:=proc(M,q,x) local n,k,gu:
k:=nops(M):
n:=nops(M[1]):

 gu:=WtEnu(M,q,x):

expand(normal((1+(q-1)*x)^n/q^k*subs(x=(1-x)/(1+(q-1)*x),gu))):

end:





#MinWtD(M,q): The minimum weight of the DUAL of a subspace of GF(q)^n generated by M (where n:=nops(M[1]). Try:
#M:=Rqnk(2,10,3); NinWtD(M,2);
MinWtD:=proc(M,q) local x:
ldegree(WtEnuD(M,q,x)-1,x):
end:


#MinWtStatD(q,n,k,N,K): inputs pos. integer q (prime or prime power), and pos. integers n and k, 1<=k<=n, and a small pos. integer N, and a large positive integer K
#uses the Calabi-Wilf algorithm to find 
#(i) upper bounds for the smallest weight of the DUAL of  a k-dim subspace of GF(q)^n
#(ii) lower bounds for the highest smallest weight in a  DUAL k-dim subspace of GF(q)^n
#(iii) estimated average, variance, and scaled moments up to the N-th. Try:
#(iv): A vector space that achives the largest minimal weight for the dual
#MinWtStatD(2,50,2,4,1000);
MinWtStatD:=proc(q,n,k,N,K) local f,x,M,i,lu,f1,si,aluf:
f:=0:
si:=1:
aluf:=0:
for i from 1 to K do
 M:=Rqnk(q,n,k):
 lu:=MinWtD(M,q):
if lu>si then
  si:=lu:
 aluf:=M:
fi:
 f:=f+x^lu:
od:

f1:=f/K:



[[ldegree(f,x),degree(f,x),evalf(Alpha(f1,x,N))],aluf]:
end:

#PeIn(M,q): The perfection index (a perfect code has value 1) for the code. If the minimal weight is even it returns FAIL. Try:
#M:=Rqnk(2,10,5); PeIn(M,2);
PeIn:=proc(M,q) local d,t,k,n,i:

d:=MinWt(M,q):
if d mod 2 =0 or d=1 then
 RETURN(FAIL):
fi:

t:=(d-1)/2:


k:=nops(M):
n:=nops(M[1]):
add(binomial(n,i)*(q-1)^i,i=1..t)/q^(n-k):
end:


#PeInD(M,q): The perfection index (a perfect code has value 1) for the code. If the minimal weight is even it returns FAIL. Try:
#M:=Rqnk(2,10,5); PeIn(M,2);
PeInD:=proc(M,q) local d,t,k,n,i:

d:=MinWtD(M,q):
if d mod 2 =0 or d=1 then
 RETURN(FAIL):
fi:

t:=(d-1)/2:


k:=nops(M):
n:=nops(M[1]):
add(binomial(n,i)*(q-1)^i,i=1..t)/q^k:
end: