### written by: Aurora Hiveley and Doron Zeilberger, Rutgers university
### last updated: July 14, 2026

print(`This is MERGING.txt, a Maple package by Aurora Hiveley and Doron Zeilberger(Rutgers University)`):
print(`accompanying their article: The Combinatorics of Lane Merging available from`):
print(``):
print(`https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/merging.html`):
print(``):
print(`Also from Aurora Hiveley's web-site:`):
print(`https://aurorahiveley.github.io/`):

print(`For a list of the functions, type Help(); for help with a specific function type: Help(FunctionName);`):

with(combinat):

Help:=proc(): 

if nargs=0 then 

print(`ALL FUNCTIONS:`):
print(`backHookCounter(T,c), CF(endstate), EVbf(k,n,i:=1), EVc(k,n,i:=1), freq(L), G2(m,n), G3(a,b,c), G4(a,b,c,d),  `): 
print(`gf(k,n,i,var:=x), gfcf(k,n,var:=x), isSYT(T), K(k,n), kLane(k,A), kVote(k,A), kVoteOld(k,A), laneLens(k,A), `): 
print(`majorityP(k,A), mergingP(k,A), multigf(k,n,var := x), numSeqToSYT(T,k), NuSYT(L), Park(n,k), SD(k,n,i:=1), `): 
print(`SYTaddn(T), SYTdelete(T), SYTdist(T), SYTgf(k,n,var := x), SYTgfbf(k,n, var := x), SYTkn(k,n), SYTshape(T), Var(k,n,i:=1)`):


## LANE SIMULATION TOOLS

elif nargs=1 and args[1]=kLane then
print(`kLane(k,A): computes the ending configuration of the lane merging scenario represented by the arrival sequence A for k lanes`);
print(`input: k-ary string representing arrival sequence A`);
print(`output: a list L = [L1,L2,..,Lk] where L1 is the list (in order) of cars in the right lane, L2 is the same for the next lane, `);
print(`ending with Lk as the same for the leftmost lane.`);
print(`example: kLane(2,[1,1,2,2,2,1,1,2]); kLane(3,[1,1,2,3,3,2,1,3,2]);`);

elif nargs=1 and args[1]=laneLens then 
print(`laneLens(k,A): the lengths of all k lanes after the arrival sequence A has been resolved`):
print(`input: arrival sequence A and number of lanes k`):
print(`output: a list of lengths of lanes after the arrival sequence has been resolved,`): 
print(`i.e., L[1] = number of cars in the rightmost lane, through L[k] = number of cars in the leftmost lane`):
print(`example: laneLens(3,[1,1,2,3,3,2,1,3,2])`):

elif nargs=1 and args[1]=K then 
print(`K(k,n): outputs the set of all strings of length n with elements {1,..,k}`);
print(`input: alphabet size k, string length n`);
print(`output: set of k-ary strings encoded as lists`):
print(`example: K(2,5); K(3,3);`):

elif nargs=1 and args[1]=kVoteOld then
print(`kVote(k,A): computes the ending configuration of the voting scenario represented by the voting sequence A for k candidates`);
print(`input: k-ary string representing voting sequence A`);
print(`output: a list L = [L1,L2,..,Lk] where L1 is the list (in order) of voters who voted for the most popular candidate, `);
print(` L2 is the same for the second most popular candidate, and ending with Lk as the same for the least popular candidate.`);
print(`example: kVote(3,[1,1,2,3,3,2,1,3,2]);  kVote(3,[1,1,1,2,2,3,3,3,3])`);

elif nargs=1 and args[1]=kVote then
print(`kVote(k,A): computes the ending configuration of the voting scenario represented by the voting sequence A for k candidates`);
print(`input: k-ary string representing voting sequence A`);
print(`output: a list L = [L1,L2,..,Lk] where L1 is the list (in order) of voters who voted for the most popular candidate, `);
print(` L2 is the same for the second most popular candidate, and ending with Lk as the same for the least popular candidate.`);
print(`example: kVote(3,[1,1,1,2,2,3,3,3,3]); `);


### GENERATING FUNCTIONS

elif nargs=1 and args[1]=gf then 
print(`gf(k,n,i:=1,var := x): the statistical generating function representing the distribution across all arrival sequences in K(k,n)`):
print(`counting the number of cars in the lane indexed i (where 1 is the rightmost lane and k is the leftmost.)`):
print(`input: number of lanes k, arrival sequence length n, lane index i (default 1), and optional parameter for variable (default x)`):
print(`output: the statistical generating function, i.e., a polynomial where the coefficient of x^m is `):
print(`the number of arrival sequences in K(k,n) which result in m cars in the lane index i`):
print(`example: gf(2,4,1); gf(3,6,1);`):

elif nargs=1 and args[1]=multigf then 
print(`multigf(k,n,var := x): the multivariate statistical generating function representing the distribution across all arrival sequences in K(k,n)`):
print(`where each term takes the form x_1^m_1*x_2^m_2*...*x_k^m_k, counting the number of arrival sequences which end with this configuration`):
print(`input: number of lanes k, arrival sequence length n, and optional parameter for variable (default x)`):
print(`output: the multivariate statistical generating function, i.e., a polynomial where the coefficient of x_1^m_1*x_2^m_2*...*x_k^m_k is `):
print(`the number of arrival sequences in K(k,n) which result with this configuration`):
print(`example: multigf(2,4); multigf(3,6);`):

elif nargs=1 and args[1]=gfCF then 
print(`gfCF(k,n,var := x): the multivariate statistical generating function representing the distribution across all arrival sequences in K(k,n)`):
print(`where each term takes the form x_1^m_1*x_2^m_2*...*x_k^m_k, counting the number of arrival sequences which end with this configuration`):
print(`input: number of lanes k, arrival sequence length n, and optional parameter for variable (default x)`):
print(`output: the multivariate statistical generating function, i.e., a polynomial where the coefficient of x_1^m_1*x_2^m_2*...*x_k^m_k is `):
print(`the number of arrival sequences in K(k,n) which result with this configuration, calculated using the closed form for all partitions of n into k pieces`):
print(`example: gfCF(2,4); gfCF(3,6);`):


### RECURRENCES AND CLOSED FORM ENUMERATION 

elif nargs = 1 and args[1]=G2 then 
print(`G2(m,n): computes the number of arrival sequences which end with m cars in the right lane and n in the left lane`): 
print(`using the recurrences from Bardenova et al. Lemma 5 `):
print(`input: integers m and n representing final lane lengths`):
print(`output number of arrival sequences whose paths end with those lane lengths`):
print(`example: G2(2,0); G2(6,5);`):

elif nargs = 1 and args[1]=G3 then 
print(`G3(a,b,c): using recurrence cases, computes the number of arrival sequences`):
print(`which end with a cars in the right lane, b in the middle, and c in the left lane`): 
print(`input: nonincreasing integers a,b,c representing the lane lengths from right to left`):
print(`output: number of arrival sequences which end with the input lane lengths`):
print(`example: G3(3,2,1); G3(6,4,2);`):

elif nargs = 1 and args[1] = freq then 
print(`freq(L): frequency notation for the list L`):
print(`input: list L`):
print(`output: list with frequency of each element, so freq(L)[i] = number of occurrences of i in L`):
print(`example: freq([1,3,2,1,3,4]);`):

elif nargs = 1 and args[1] = CF then 
print(`CF(endstate): calculates the number of arrival sequences which end with the configuration endstate`):
print(`using the closed form from Theorem 13`):
print(`input: a list of weakly decreasing lane lengths, endstate, where nops(endstate) = k is the number of lanes`):
print(`output: the number of arrival sequences ending with the endstate configuration, `):
print(`example: CF([2,1]); CF([4,3,2]);`):


### STANDARD YOUNG TABLEAU 

elif nargs = 1 and args[1]=SYT then 
print(`SYT(T): the set of arrival sequences of length n (= maximum element of T) with k lanes (k is the number of rows of T)`):
print(`whose ending configuration is precisely the SYT, T`):
print(`input: standard young tableau, T`):
print(`output: set of arrival sequences which result in the tableau T`):
print(`example: SYT([[1,2],[3,4],[5,6]]); SYT([[1,2,3],[4,5,6]])`):

elif nargs = 1 and args[1]=NuSYT then 
print(`NuSYT(L): the number of SYT of shape L`):
print(`input: standard young tableau shape, L, as a list of weakly decreasing positive integers`):
print(`output: number of SYT of shape L`):
print(`example: NuSYT([3,3,3]); NuSYT([2,1])`):

elif nargs = 1 and args[1] = isSYT then 
print(`isSYT(T): returns boolean for if the shape T is a SYT. checks legality of shape, increasing rows and columns, and standard-ness`):
print(`input: shape T, given as a list of lists corresponding to each row of the tableau`):
print(`output: boolean true/false `):
print(`example: isSYT([[1,2,3],[4,5],[6]]); isSYT([[1,2],[3,4,5],[6]]); isSYT([[1,2,4],[4,5],[6]])`):

elif nargs = 1 and args[1]=SYTdist then 
print(`SYTdist(T): the list of possible entries, at each index, of an arrival sequence which produces the tableau T`):
print(`input: standard young tableau, T`):
print(`output: list of sets corresponding to the allowed entries of the sequence at that position`):
print(`example: SYTdist([[1,2],[3,4],[5,6]]); SYTdist([[1,2,3],[4,5,6]]);`):

elif nargs = 1 and args[1]=SYTdelete then 
print(`SYTdelete(T): recursive SYT resulting from deletion of largest element n, and the number of possible candidate`):
print(`digits in {0,..,k-1} which could have been in position "n" of an arrival sequence which produced the tableau T`):
print(`input: standard young tableau, T`):
print(`output: new SYT after deletion of cell labelled n, and number of possible digits in entry n of arrival sequence`):
print(`example: SYTdelete([[1,2],[3,4],[5,6]]); SYTd([[1,2,6],[3,4],[5]]);`):

elif nargs = 1 and args[1]=SYTaddn then 
print(`SYTaddn(T): set of SYT acquired by adding next largest element to the existing tableau T in each legal position`):
print(`IMPORTANTLY: while maintaining the number of (possibly empty) rows of the tableau`):
print(`input: standard young tableau, T`):
print(`output: set of new SYT obtained by adding n in each legal position, maintaining existing number of rows`):
print(`example: SYTaddn([[1,2],[3,4],[5]]); SYTaddn([[1,2],[3],[]]);`):

elif nargs = 1 and args[1] = SYTkn then 
print(`SYTkn(k,n): set of all SYT with n entries and at most k rows`):
print(`input: integers k and n `):
print(`output: set of all SYT with n entries and k rows (possibly empty)`):
print(`example: SYTkn(2,4); SYTkn(3,5)`):

elif nargs = 1 and args[1] = SYTshape then 
print(`SYTshape(T): returns the shape of the SYT T as a list of (weakly decreasing) row lengths`):
print(`input: tableau T `):
print(`output: list of row lengths from top to bottom`):
print(`example: SYTshape([[1,2,3,4],[5,6],[7]]);`):


### SYT MAPPING ENUMERATION 

elif nargs = 1 and args[1] = backHookCounter then 
print(`backHookCounter(T,c): the number of cells in the column to the left of cell c which are below cell c`):
print(`and which are smaller than the entry in cell c`):
print(`input: a SYT, T, and the location of a cell c = [i,j]`):
print(`output: integer number of cells in the backwards oriented hook from c which are smaller than the entry in cell c`):
print(`example: backHookCounter([[1,4,6],[3,5],[7]], [1,2]);  backHookCounter([[1,4,6],[3,5],[7]], [2,2]); `)

elif nargs = 1 and args[1] = SYTgfbf then 
print(`SYTgfbf(k,n,var := x): statistical generating function of arrival sequences as SYT with k rows and n elements `):
print(`input: number of lanes k and number of tableau entries n, optional parameter for variable var (default x)`):
print(`output: the statistical generating function, computed by brute force, where the coefficient `):
print(`of var[T] for a tableau T is the number of arrival sequences whose final arrangement is the tableau T `):
print(`example: SYTgfbf(2,4); `): 

elif nargs = 1 and args[1] = SYTgf then 
print(`SYTgfbf(k,n,var := x): statistical generating function of arrival sequences as SYT with k rows and n elements `):
print(`using the product of back hook lengths`):
print(`input: number of lanes k and number of tableau entries n, optional parameter for variable var (default x)`):
print(`output: the statistical generating function, computed by the conjectured product formula, where the coefficient `):
print(`of var[T] for a tableau T is the number of arrival sequences whose final arrangement is the tableau T `):
print(`example: SYTgfbf(2,4); `): 

elif nargs = 1 and args[1] = numSeqToSYT then 
print(`numSeqToSYT(T,k): value computed by conjectured product formula for the number of arrival sequences`):
print(`which resolve to the tableau T`):
print(`input: tableau T and number of lanes k `):
print(`output: the number of arrival sequences for k lanes which produce the final arrangement in `):
print(`the tableau T, computed with the conjectured product formula`):
print(`example: numSeqToSYT([[1,3,5],[2,4]],2)`):

elif nargs = 1 and args[1] = ff then 
print(`ff(n,k): falling factorial n^\underline{k} `):
print(`input: integers n and k (k <= n)`):
print(`output: falling factorial n^\underline{k} = n*(n-1)*...*(n-k+1)`):
print(`example: ff(4,2);`):


### WALKS 

elif nargs=1 and args[1]=mergingP then 
print(`mergingP(k,A): generates the merging path of points p produced by a k-ary arrival sequence A`); 
print(`according to a generalization of the framework outlined in Bardenova et al. (2022)`):
print(`now, coordinate p1 is incremented when a car joins the rightmost lane, coordinate p2 is incremented`):
print(`when a car joins the next lane, and so on. finally, coordinate pk is incremeneted when a car joins the rightmost lane.`):
print(`input: number of lanes k and a k-ary arrival sequence A`):
print(`output: a list of points [p1,p2,...,pk] hit by the merging path`):
print(`example: mergingP(2,[1,1,2,2,2,1,1,2]); mergingP(3,[1,1,2,3,3,2,1,3,2]);`):

elif nargs=1 and args[1]=majorityP then 
print(`majority(k,A): generates the majority path of points p produced by a voting sequence with k choices`); 
print(`where a point p = [p1, p2, ..., pk] is the number of votes of the most to least popular choices`):
print(`so p[i] is the number of votes for the i-th most popular choice after m := add(p[i], i=1..k) votes have been cast`);
print(`input: number of choices k and a voting sequence A`):
print(`output: a list of points [p1,p2,...,pk] hit by the majority path`):
print(`example: majorityP(2,[1,1,2,2,2,1,1,2]); majorityP(3,[1,1,2,3,3,2,1,3,2]);`):



### STATISTICS 

elif nargs=1 and args[1]=EV then 
print(`EV(k,n,i:=1): computes the expected length of the lane indexed by i for k lanes and arrival sequences of l`):
print(`input: number of lanes k, lane index i (from 1 to k, counted right to left, default is 1), and arrival sequence length l`):
print(`output: expected length of lane i. multiply by k^l to count generalization of R(B_l)`):
print(`example: EV(2,1,4); EV(3,3,5);`):

elif nargs=1 and args[1]=EVc then 
print(`EVc(k,n,i:=1): computes the expected length of the lane indexed by i for k lanes and arrival sequences of l`):
print(`input: input the number of lanes k, number of cans n, and index of lane of interest i (from right to left, default 1)`):
print(`output: expected length of lane i, calculated cleverly using the closed form`):
print(`example: EVc(2,1,4); EVc(3,3,5);`):

elif nargs=1 and args[1]=Varbf then 
print(`Var(k,n,i:=1): calculates, using the proven closed form, the variance of the length of lane i over all k-lane arrival sequences of length n`):
print(`input: number of lanes k, number of cars n, and the index i of lane of interest (1 <= i <= k, default is 1)`):
print(`output: variance, as a decimal`):
print(`example: Var(4,2,1);`):

elif nargs=1 and args[1]=SD then 
print(`SD(k,n,i:=1): calculates, using the proven closed form, the standard deviation of the length of lane i over all k-lane arrival sequences of length n`):
print(`input: number of lanes k, number of cars n, and the index i of lane of interest (1 <= i <= k, default is 1)`):
print(`output: standard deviation, as a decimal`):
print(`example: SD(4,2,1);`):



else 
print(`No such function exists. Type Help(); for a list of functions.`):

fi:

end:




## K LANE SIMULATION 
# 
# input: number of lanes k and a k-ary string A corresponding to an arrival sequence
# output: list L = [L1,...,Lk] where Li is the list (in order) of cars in lane i 
# note that the lanes are read from right to left, so L[1] is the rightmost lane, and L[k] is the leftmost  

kLane := proc(k,A) local i,j,m,L,cons,consR:
if not convert(A,set) subset {seq(1..k)} then 
    error(`Car preferences require more than the input k many lanes.`):
fi:

L := [seq([], i=1..k)]:     # initialize lane lists 

for j from 1 to nops(A) do 
    cons := L[(1..A[j])]:  # lanes to consider based on car preference 
    m := min[index]([seq(nops(cons[i]), i=1..nops(cons))]):     # index of shortest lane, with the first lane (rightmost) selected in a tie 
    L := [ op(L[1..m-1]), [op(L[m]),j], op(L[m+1..-1])]:
od:

L: # note that output is right lane to left lane 
end:


## K CANDIDATE VOTING SIMULATION (OLD)
# 
# input: number of candidates k and a k-ary string A corresponding to the votes cast
# output: list L = [L1, ..., Lk] where Li is the list (in order) of individuals who votes for the i-th most popular candidate 

kVoteOld := proc(k,A) local i,j,V,v:
V := []:
v := {}:

for i from 1 to k do                # fix a candidate
    for j from 1 to nops(A) do      # check all voters
        if A[j] = i then 
            v := v union {j}:       # voters who voted for candidate i 
        fi:
    od: 

    V := [op(V), convert(v,list)]:
    v := {}:                        # reset 
od:

V := sort(V, `=`(key, x -> nops(x))):  # sort by row lengths, i.e. by candidate popularity 
[seq(V[nops(V)-i+1], i=1..nops(V))]:   # reverse order so takes form of YT (weakly decreasing row lengths)
end:



## K CANDIDATE VOTING SIMULATION
# 
# input: number of candidates k and a k-ary string A corresponding to the votes cast
# output: tableau formed by adding element n to row i when voter n casts a vote for the i-th most popular candidate  

kVote := proc(k,A) local i,j,T,votes,numVotes,sortedNumVotes,c,r:
T := [[]$k]: # initialize tableau
votes := [[]$k]:

for i from 1 to nops(A) do 
    c := A[i]:
    votes := [op(votes[1..c-1]), [op(votes[c]), i], op(votes[c+1..-1])];
    numVotes := [seq(nops(votes[i]), i=1..k)]:
    sortedNumVotes := sort(numVotes):

    # which rows of the SYT have the number of votes that candidate c has 
    # put i in the lowest row (first occurring) with that many votes?
    member(nops(votes[c]), sortedNumVotes, 'r'):
    T := [op(T[1..r-1]), [op(T[r]),i], op(T[r+1..-1])]:
od:

[seq(T[k-i], i=0..k-1)]:   # reverse the order 
end:


## K-ARY STRINGS 
# 
# input string length n and size of alphabet k 
# output set of all strings of length n with entries 1,...,k

K := proc(k,n) local K1,K2,kappa,k1:
option remember:

if n = 1 then 
    RETURN({ seq([kappa], kappa=1..k)}):
else 
    K1 := K(k,n-1):
    K2 := { seq( seq( [op(k1),kappa] , kappa=1..k ), k1 in K1) }:
fi: 

K2:
end:


## LANE LENGTHS 
# 
# input arrival sequence A and number of lanes k 
# output a list of lengths of lanes after the arrival sequence has been resolved, i.e., L[1] = number of cars in the rightmost lane, 
# through L[k] = number of cars in the leftmost lane 

laneLens := proc(k,A) local i,L:
L := kLane(k,A):
[seq(nops(L[i]), i=1..k)]:
end:


## STATISTICAL GENERATING FUNCTION FOR K LANES 
# 
# input number of lanes k, arrival sequence length l, lane index i, and optional parameter for variable (default x)
# output polynomial where the coefficient of x^n is the number of arrival sequences which result in n cars in the lane index i 

gf := proc(k,n,i:=1,var := x) local A:
if not i >= 1 and i <= k then 
    error(`Lane index not in range {1,..,k}`):
fi:

add(var^nops(kLane(k,A)[i]), A in K(k,n)):
end:


## MULTIVARIATE STATISTICAL GENERATING FUNCTION FOR K LANES 
# 
# input number of lanes k, arrival sequence length l, and optional parameter for variable (default x)
# output polynomial where the coefficient of x_1^n_1*x_2^n_2*...*x_k^n_k is the number of arrival sequences which result in [n_1, n_2, ..., n_k] 

multigf := proc(k,n,var := x) local A,a,f,i:
# to cut down on runtime, but does increase memory 
f := 0: 

for A in K(k,n) do 
    a := kLane(k,A):
    f := f + mul( var[i]^nops(a[i]) , i=1..k ):
od: 

f:
end:


## GENERATING FUNCTION USING CLOSED FORM 
# 
# input number of lanes k, number of cars n, and optional variable (default x)
# outputs the multivariate generating function calculated using closed form for all partitions of n into k pieces

gfCF := proc(k,n,var:=x) local p,i:
add(CF(p)*mul(var[i]^p[i],i=1..k) , p in Parnk(n,k)):
end:



## EXPECTED VALUE BY BRUTE FORCE 
# 
# input number of lanes k, number of cars n, and index of lane of interest i (from right to left, default 1)
# output expected number of cars in lane indexed by i for a fixed length l 

EVbf := proc(k,n,i:=1) local A:
add(nops(kLane(k,A)[i]), A in K(k,n))/k^n:
end:



### EXPECTED VALUE COMPUTED CLEVERLY 
# 
# input the number of lanes k, number of cans n, and index of lane of interest i (from right to left, default 1)
# output expected number of cars in lane indexed by i, using closed form from Theorem 13

EVc := proc(k,n,i:=1) local f,g,j:
f := gfCF(k,n,x):
g := diff(f,x[i]):
subs({seq(x[j]=1, j=1..k)}, g)/(k^n):
end:



## 2 LANE RECURRENCE 
# 
# input integers m >= n representing number of cars in ending state with n cars in left lane, m cars in right lane 
# output number of arrival sequences which end with those lane lengths, using recurrence 

G2 := proc(m,n)
option remember: 

if m = 0 and n = 0 then 
    RETURN(1):
elif m <> 0 and n = 0 then 
    RETURN(2):
fi:

if m > n+1 then 
    G2(m,n-1) + G2(m-1,n):
elif m = n+1 then 
    G2(m,n-1) + 2*G2(m-1,n):
elif m = n then
    G2(m,n-1):
else 
    RETURN(0):   # m < n --> not allowed!
fi:

end:


### G3 USING RECURRENCE 
# 
# input positive integers a,b,c 
# output the number of arrival sequences which resolve to the SYT with shape [a,b,c] using recurrence  

G3 := proc(a,b,c) 
option remember: 

if not (a >= b and b >=c and c >= 0) then 
    RETURN(0):
elif a=0 and b=0 and c=0 then 
    RETURN(1):
fi: 

if a-b > 1 and b-c> 1 then 
    G3(a-1,b,c) + G3(a,b-1,c) + G3(a,b,c-1):     # most common case 
elif a-b <= 1 and b-c > 1 then 
    2*G3(a-1,b,c) + G3(a,b-1,c) + G3(a,b,c-1):   # case 2 & 3 together 
elif a-b > 1 and b-c <= 1 then 
    G3(a-1,b,c) + 2*G3(a,b-1,c) + G3(a,b,c-1):   # case 4 & 5 together
elif a-b =1 and b-c = 1 then 
    2*G3(a-1,b,c) + 2*G3(a,b-1,c) + G3(a,b,c-1):
else                                                # at least one of a-b or b-c is = 0 
    3*G3(a-1,b,c) + 2*G3(a,b-1,c) + G3(a,b,c-1):
fi:

end:


### G4 RECURRENCE 
# 
# input positive integers a,b,c,d
# output the number of arrival sequences which resolve to the SYT with shape [a,b,c] using recurrence

G4 := proc(a,b,c,d) local f,i,L,L1,D,t,co:
option remember: 

if not (a >= b and b >=c and c >= d and d >= 0) then           
    RETURN(0):
elif a=0 and b=0 and c=0 and d=0 then 
    RETURN(1):
fi: 

f := 0:     # recurrence function, will modify as we go 
L := [a,b,c,d]: 
D := [seq(L[i]-L[i+1], i=1..3)]:

for i from 1 to 3 do 
    L1 := L - [0$(i-1),1,0$(4-i)]:         
        
    if L[i] - L[i+1] > 1 then  
        f := f + G4(op(L1)):

    elif L[i] - L[i+1] = 0 then 
        f := f : 
        
    elif L[i] - L[i+1] = 1 then 
        co := 2:

        for t from i+1 to 4-1 while L[t] - L[t+1] = 0 do 
            co ++:
        od: 

        f := f + co*G4(op(L1)):

    fi:
od:

f := f + G4(op(L - [0$(4-1),1])):
 
end:



### PARTITIONS OF N INTO K (POSSIBLY EMPTY) PIECES
#
# input integers n and k 
# outputs the set of weakly decreasing nonnegative sequences of length k that add up to k (allowing zero!)

Parnk:=proc(n,k) local a:
option remember:
{seq(op(Parnk1(n,k,a)),a=0..n)}:
end:


## HELPER: RECURSIVE CONSTRUCTOR OF PARTITIONS
# 
# input integers n, k, and a 
# outputs the set of weakly decreasing nonnegative sequences of length k that add up to k with largest part a

Parnk1:=proc(n,k,a) local S,a1,s,S1:
option remember:
if k=1 then
if n=a then
 RETURN({[a]}):
 else
  RETURN({}):
fi:
fi:

S:={}:

for a1 from 0 to a do
 S1:=Parnk1(n-a,k-1,a1):
 S:=S union {seq([a, op(s)], s in S1)}:
od:
S:
end:


### CLOSED FORM 
# 
# input a list of lane lengths corresponding to an endstate 
# outputs the number of k-lane arrival sequences ending in that configuration by the closed form from Theorem 13

CF := proc(endstate) local p,i,e,k:
k := nops(endstate): 

for i from 1 to k-1 do 
    if not endstate[i] >= endstate[i+1] then 
        RETURN(0):
    fi:
od:

p := Afreq(endstate):
multinomial(k,op(p))*multinomial(add(e, e in endstate), op(endstate))

end:


### SYT OF SHAPE L 
# copied from code written in 642:640 at rutgers university in spring 2026 
# 
# input a SYT shape L as a list of nonincreasing positive integers 
# output the set of all SYT with the input shape 

SYT:=proc(L) local i, k,n,S,L1,S1,s1,s,LL,kk,ext:
option remember:

k:=nops(L):
n:=add(L[i],i=1..k):
if k=0 then
    RETURN({[]}):
fi:

# handle empty rows 
if member(0,L) then 
    member(0,L,'i'):
    LL := L[1..i-1]: # slice off the 0s 
    ext := k - (i-1): # number of empty rows to add back on at end 
else 
    LL := L:
    ext := 0:
fi:

kk := k-ext: # number of rows being used

S:={}:
# search for all rows where entry n can be legally placed 
for i from 1 to kk-1 do
    if LL[i]>LL[i+1] then
        L1:=[op(1..i-1,LL),LL[i]-1,op(i+1..kk,LL)]:
        S1:=SYT(L1):
        S:=S union {seq([op(1..i-1,s1),[op(s1[i]),n],op(i+1..kk,s1)],s1  in S1)}:
    fi:
od:

if LL[kk]>1 then
    L1:=[op(1..kk-1,LL),LL[kk]-1]:

    S1:=SYT(L1):
    S:=S union {seq([op(1..kk-1,s1),[op(s1[kk]),n]],s1  in S1)}:
else
    L1:=[op(1..kk-1,LL)]:
    S1:=SYT(L1):
    S:=S union {seq([op(1..kk-1,s1), [n]],s1  in S1)}:
fi:

# add zeros back on
{seq( [op(s),[]$ext], s in S)}:
end:


### NUMBER OF SYT OF SHAPE L 
# copied from code written in 642:640 at rutgers university in spring 2026 
# 
# input a SYT shape L as a list of nonincreasing positive integers 
# output the number of all SYT with the input shape (i.e., nops(SYT(L)))

NuSYT:=proc(L) local i, k,S,L1:
option remember:
k:=nops(L):

if k=0 then
    RETURN(1):
fi:

S:=0:

# search for all rows where entry n can be legally placed 
for i from 1 to k-1 do
    if L[i]>L[i+1] then
        L1:=[op(1..i-1,L),L[i]-1,op(i+1..k,L)]:
        S:=S+NuSYT(L1): 
    fi:
od:

if L[k]>1 then
    L1:=[op(1..k-1,L),L[k]-1]:
    S:=S+ NuSYT(L1):

else
    L1:=[op(1..k-1,L)]:
    S:=S+NuSYT(L1):
fi:

S:
end:


## ALL SYT OF SIZE n WITH k ROWS 
# 
# input integers k and n
# outputs the set (generated recursively) of all SYT with k rows and n entries

SYTkn := proc(k,n) local T: 
option remember: 

if n=0 then 
    RETURN ({[[]$k]}):
fi:

{seq( op(SYTaddn(T)) , T in SYTkn(k,n-1))}: 
end:


## HELPER: SYT ADD ELEMENT "n" IN ALL POSSIBLE POSITIONS 
#
# inputs current tableau T 
# outputs the set of all tableau (WITH SAME NUMBER OF ROWS AS T) with element n added to each legal position

SYTaddn := proc(T) local n,k,S,i:

if convert(T,set) = {[]} then 
    n := 1:
else
    n := max(T) +1:    # next largest element to add to SYT
fi: 

k := nops(T):

S := { [ [op(T[1]),n], op(T[2..-1]) ] }:  # add n to the end of the first row (always allowed)

# add element n to the end of any existing row, or to its own new row! full shape cannot have more than k rows 
for i from 2 to k do 
    if nops(T[i])+1 <= nops(T[i-1]) then   
        S := S union { [ op(T[1..i-1]), [op(T[i]),n], op(T[i+1..k])  ] }:
    fi:
od:

S:
end:


## SYT DISTRIBUTION 
# 
# input a standard young tableau T 
# output the list of choices for each arrival sequence entry in an arrival sequence which produces the SYT, T 

SYTdist := proc(T) local S,n,L,i,j:
if not isSYT(T) then 
    error `Input tableau is not a SYT`: 
fi:

S := convert(SYT(T),list):
n := max(T):
L := []:

for i from 1 to n do 
    L := [op(L), {seq(S[j,i], j=1..nops(S))}]:
od:

L:
end:


## SYT RECURSIVE DELETION OF BOX "n"
#
# input a SYT which represents a final configuration for a lane merging problem 
# outputs a new SYT, T1, which removes the box labelled n, and the number of possible candidate digits in {0,..,k-1}
# which could have been in position "n" of an arrival sequence which produced the tableau T

SYTdelete := proc(T) local n,k,c,r,T1,count: 
if not isSYT(T) then 
    error `Input tableau is not a SYT`: 
fi:

n := max(T):
k := nops(T): 
c := [max[index](T)]:     # location of the maximum element, n

# new tableau with element n removed 
T1 := [op(T[1..c[1]-1]), [op(T[c[1]][1..-2])], op(T[c[1]+1..-1])]:

# how many open spots are below the position c when the number n was placed? equivalent to how many digit options we have in the arrival sequence
count := 1:

for r from c[1]+1 to k while nops(T[r]) = c[2]-1 do
    count++:
od:

count, T1:
end:



### IS AN SYT (ERROR CHECKER )
# 
# inputs a tableau T 
# outputs boolean true/false to the question "is T a standard young tableau?" 

isSYT := proc(T) local n,k,i,j,lens:
k := nops(T):
n := add(nops(T[i]), i=1..k):

# check shape 
lens := [seq(nops(T[i]), i=1..k)]:
for j from 2 to nops(lens) do 
    if not lens[j-1] >= lens[j] then 
        RETURN(false):
    fi:
od:

# check *standard* young tableau
if not {seq( op(T[i]) , i=1..k )} = {seq(1..n)} then      # set of entries is not {1,...,n = number of boxes}
    RETURN(false):
fi:

# check increasing rows 
for i from 1 to k do 
    for j from 2 to nops(T[i]) do 
        if not T[i][j-1] < T[i][j] then
            RETURN(false):
        fi:
    od:
od:


# check increasing columns. j is the column number
for j from 1 to nops(T[1]) do 
    for i from 2 to k while nops(T[i]) >= j do      # for each successive row whose length is long enough to have an entry
        if not T[i-1][j] < T[i][j] then 
            RETURN(false):
        fi:
    od:
od:

true:
end:

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

### SYT MAPPING ENUMERATION 

## NUMBER OF SEQUENCES MAPPING TO AN SYT
# 
# input tableau T and number of lanes k 
# output the number of arrival sequences for k lanes which produce the final arrangement in 
# the tableau T, computed with the product formula 

numSeqToSYT := proc(T,k) local i,j,d,m:
if not isSYT(T) then 
    error `Input tableau is not a SYT`: 
fi:

if not nops(T) <= k then 
    error `Tableau has too many rows for suggested number of lanes.`
fi:

# number of occupied rows, m 
m := nops( convert(T,set) minus {[]} ):   # every nonempty row should be distinct, since *is* a SYT! 

ff(k,m)*mul( mul(1+backHookCounter(T,[i,j]), j=2..nops(T[i]) ), i=1..k-1 ):
end:


## NUMBER OF SEQUENCES MAPPING TO AN SYT - CLEVER
# same as conjprod, but uses recursion adapted from SYTdelete to speed up computation time 
# 
# input tableau T and number of lanes k 
# output the number of arrival sequences for k lanes which produce the final arrangement in 
# the tableau T, computed with the product formula 

numSeqToSYTc := proc(T,k) local i,j,d:
option remember: 

if not isSYT(T) then 
    error `Input tableau is not a SYT`: 
fi:

if not nops(T) <= k then 
    error `Tableau has too many rows for suggested number of lanes.`
fi:

sts1(T):
end:


## HELPER: RECURSIVE PRODUCT CALCULATOR/TABLEAU DELETION
#
# input tableau T 
# output product obtained from multiplying the number of possible digits in a k-ary arrival sequence 
# which resolves to the SYT T

sts1 := proc(T) local n,k,c,r,T1,count:
option remember: #?

if convert(T,set) = {[]} then 
    RETURN(1):
fi:

n := max(T):
k := nops(T): 
c := [max[index](T)]:     # location of the maximum element, n

# new tableau with element n removed 
T1 := [op(T[1..c[1]-1]), [op(T[c[1]][1..-2])], op(T[c[1]+1..-1])]:

# how many open spots are below the position c when the number n was placed? equivalent to how many digit options we have in the arrival sequence
count := 1:

for r from c[1]+1 to k while nops(T[r]) = c[2]-1 do
    count++:
od:

count*sts1(T1):
end:


## BACKWARDS HOOK COUNTER 
# 
# input tableau T and location of a cell c (as coordinates [i,j])
# output the number of cells in the column j-1 and row i+1 or larger which are smaller than the entry in cell c

backHookCounter := proc(T,c) local i,j,k,x,A,count:
i := c[1]:
j := c[2]:
k := nops(T):
count := 0:

A := T[i,j]:    # element at top of hook being considered

for x from i+1 to k while nops(T[x]) >= j-1 do 
    if A > T[x][j-1] then 
        count++
    fi:
od:

count:
end:


## FALLING FACTORIAL 
#
# input integers n and k 
# output falling factorial n^\underline{k}

ff := proc(n,k) local i:
mul(n-i, i=0..k-1):
end:

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

## SYT GENERATING FUNCTION - BY BRUTE FORCE 
#
# input number of lanes k and number of tableau entries n 
# output the statistical generating function, computed by brute force, where the coefficient 
# of var[T] for a tableau T is the number of arrival sequences whose final arrangement is the tableau T 

SYTgfbf := proc(k,n, var := x) local A:
add(var[kLane(k,A)],A in K(k,n)) ;
end:


## SYT SHAPE 
# 
# input a SYT, T 
# output its shape (list of row lengths from top to bottom)

SYTshape := proc(T) local i: 
[seq(nops(T[i]), i=1..nops(T))]:
end:


## SYT GENERATING FUNCTION - USING PRODUCT
#
# input number of lanes k and number of tableau entries n 
# output the statistical generating function, computed from our back hook product formula, 
# where the coefficient of var[T] for a tableau T is the number of arrival sequences whose 
# final arrangement is the tableau T 

SYTgf := proc(k,n,var := x) local T:
add( conjprod(T,k)*var[T] , T in SYTkn(k,n)):
end:


## FREQUENCY VECTOR 
#
# inputs a list L, outputs a list containing frequency of each l in L 
# so F[i] = frequency of i in L

freq := proc(L) local l,F: 
F := [0$max(L)]:

for l in L do 
    F := [op(1..l-1,F), F[l]+1, op(l+1..-1,F)]:
od:

F: 
end:


## PARTITIONS OF n INTO k PIECES 
# 
# input integers n and k 
# output set of all partitions of n into exactly k (possibly empty) parts
# i.e., the set of SYT shapes with k rows and n elements 

Park:=proc(n,k) local S,k1,S1,s1:
option remember:

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

if k=0 then
 if n=0 then
   RETURN({[]}):
   else
    RETURN({}):
fi:
fi:

S:={}:
for k1 from 0 to k do
 S1:=Park(n-k,k1):
 S:= S union {seq([op(s1+[1$k1]),1$(k-k1)], s1 in S1)}:
od:
S:
end:


## MERGING PATH FOR K LANES 
# 
# input number of lanes k and arrival sequence A 
# outputs the list of points in Z^k traced by the merging path drawn by A 
# using the framework from bardenova et al (2022)
# i.e., when a car is added to lane n, the coordinate of p corresponding to that lane is incremented by one 
# where p1 is rightmost lane, p2 is second lane, etc. 
# note that the lanes are ordered in a mirrored fashion from the ordering used in the lane merging simulations 

mergingP := proc(k,A) local p,q,P,i,j,m: 
if not convert(A,set) subset {seq(1..k)} then 
    error(`Car preferences require more than the input k many lanes.`):
fi:

p := [0$k]:
P:= [p]:

for j from 1 to nops(A) do 
    m := min[index](p[1..(A[j])]):      # shortest lane from those considered by the car, with the first lane (rightmost) selected in a tie 
    p := [op(p[1..m-1]), p[m]+1, op(p[m+1..-1])]:
    P := [op(P),p]:
od:

P:
end:


## MAJORITY PATH FOR K LANES 
# 
# input number of choices k and voting sequence A 
# outputs the list of points in Z^k traced by the majority path drawn by A 
# a point p = [p1, p2, ..., pk] indicates that after m := add(p[i], i=1..k) votes have been cast, 
# there are p1 votes for the most popular choice, p2 for the second most popular, and so on 

majorityP := proc(k,A) local p,P,j,v: 
if not convert(A,set) subset {seq(1..k)} then 
    error(`Voter preferences require more than the input k many options.`):
fi:

p := [0$k]:     # votes for choices, in order
P:= [p]:

for j from 1 to nops(A) do 
    v := A[j]:                                # vote from person j, fixes indexing so 0 occurs first 
    p := [op(p[1..v-1]), p[v]+1, op(p[v+1..-1])]: 
    P := [op(P),sort(p, `>`)]:                      
od:

P:
end:



#### STATISTICS 

# order statistics, written by dr z during meeting 4/27 
OS := proc(v) local i,f:
f := add(x[v[i]], i=1..nops(v)):
sort( [ seq(coeff(f,x[i],1), i=0..max(v) )], `>`):
end:

# variance 
Var := proc(k,n,i:=1) local f,g,h,av,j: 
f := clevergf(k,n,x):
g := diff(f,x[i]):
av := subs({seq(x[j]=1, j=1..k)}, g)/(k^n):

h := subs({seq(x[j]=1, j=1..i-1), seq(x[j]=1, j=i+1..k)}, f):
evalf( add( coeff(h,x[i],j)*(j-av)^2, j=1..n)/(k^n) ):
end:


# standard deviation 
SD := proc(k,n,i:=1) sqrt(Var(k,n,i)): end:


#start New stuff by Dr. Z.

New:=proc(): print(` LaneSeqs(T), LaneSeqsC(T),      MulTab(T), AllSYT(k,n), RatKc(T),   VoteSeqsC(T)     `):end:

New1:=proc(): print(`GGFv(k,n,X), GGFl(k,n,x),    AllSYT(k,n), KidT(T), RatK(T), VoteTable(k,n), VoteSeqs(T), SV(v,k), Rank11(L,S), Rank1(L), ExtVo(k,v,g), LaneRange(k,v) , VoteRange(k,v) `):end:

GGFv:=proc(k,n,X) local S,s:
option remember:
S:=K(k,n):
add(X[kVote(k,s)], s in S):
end:

GGFl:=proc(k,n,X) local S,s:
option remember:
S:=K(k,n):
add(X[kLane(k,s)], s in S):
end:


#MulTab(T): Given a Young tableau finds the number of voting sequences (also lane-preferance) that lead to it
MulTab:=proc(T) local i, k,n,X:
k:=nops(T):
n:=add(nops(T[i]),i=1..k):
coeff(GGFv(k,n,X),X[T],1):
end:

AllSYT:=proc(k,n) local S,s:
S:=K(k,n):
{seq(kVote(k,s), s in S)}:
end:

#KidT(T): The child of the Young Tableau T
KidT:=proc(T) local n,k,i:
k:=nops(T):
n:=add(nops(T[i]),i=1..k):

for i from 1 to k do
if nops(T[i])<>0 and T[i][-1]=n then
RETURN([op(1..i-1,T),  [op(1..nops(T[i])-1,T[i])],op(i+1..k,T)]):
fi:
od:

end:


RatK:=proc(T):MulTab(T)/MulTab(KidT(T)): end:


#VoteTable(k,n): A table that for each tableau with n boxes and k rows (some possibly empty) that for each tableau shows
#the set of voter sequences that lead to that tableau. Try:
#VoteTable(3,5);
VoteTable:=proc(k,n) local S,T,s,V,v:
option remember:
S:=AllSYT(k,n):

for s in S do
T[s]:={}:
od:

V:=K(k,n) :

for v in V do
T[kVote(k,v)]:=T[kVote(k,v)] union {v}:
od:

[op(T),S]:
end:

#LaneTable(k,n): A table that for each tableau with n boxes and k rows (some possibly empty) that for each tableau shows
#the set of parking-preferance sequences that lead to that tableau. Try:
#LaneTable(3,5);
LaneTable:=proc(k,n) local S,T,s,V,v:
option remember:
S:=AllSYT(k,n):

for s in S do
T[s]:={}:
od:

V:=K(k,n) :

for v in V do
T[kLane(k,v)]:=T[kLane(k,v)] union {v}:
od:

[op(T),S]:
end:



#VoteSeqs(T): All the voting sequences that lead to the tableau T. Try:
#VoteSeqs([[1,2],[3,4]]);
VoteSeqs:=proc(T) local k,n,lu,i:
k:=nops(T):
n:=add(nops(T[i]),i=1..k):
lu:=VoteTable(k,n):

if not member(T,lu[2]) then
 RETURN(FAIL):
else
RETURN(lu[1][T]):
fi:
end:

#LaneSeqs(T): All the voting sequences that lead to the tableau T. Try:
#LaneSeqs([[1,2],[3,4]]);
LaneSeqs:=proc(T) local k,n,lu,i:
k:=nops(T):
n:=add(nops(T[i]),i=1..k):
lu:=LaneTable(k,n):

if not member(T,lu[2]) then
 RETURN(FAIL):
else
RETURN(lu[1][T]):
fi:
end:





#SV(v,k): The score of each party (1 through k) in the vector v
SV:=proc(v,k) local T,i:
for i from 1 to k do
T[i]:=0:
od:

for i from 1 to nops(v) do
T[v[i]]:=T[v[i]]+1:
od:
[seq(T[i],i=1..k)]:
end:


#Rank11(L,S): Given a list of integers L, and a subset of places finds the set of members of S that have the largest. Try:
#Rank11([2,3,3],{1,2,3});
Rank11:=proc(L,S) local m,S1,s:
m:=max(seq(L[s], s in S)):

S1:={}:

for s in S do
 if L[s]=m then
  S1:=S1 union {s}:
 fi:
 od:

[S1,S minus S1]:

end:

Rank1:=proc(L) local S,R,i,gu:
R:=[]:

S:={seq(i,i=1..nops(L))}:

while S<>{} do
gu:=Rank11(L,S):
R:=[op(R),gu[1]]:
S:=gu[2]:
od:
R:
end:




#ExtVo(k,v,g): Given a pos. integer k, a word v in {1, ...,k} and g between 1 and k, finds all extensions in the g-th most popular party
ExtVo:=proc(k,v,g) local L,i,S,i1,s:
L:=Rank1(SV(v,k)):
 for i from 1 to nops(L) while add(nops(L[i1]),i1=1..i)<g do od:
 S:=L[i]:
{seq([op(v),s],s in S)}:
end:


#VoteSeqsC(T): a clever (recursive) way of computing VoterSeqs(T).
VoteSeqsC:=proc(T) local k,i,n,T1,gb,V,v:
option remember:
k:=nops(T):

if T=[[]$k] then
 RETURN({[]}):
fi:
n:=add(nops(T[i]),i=1..k):

#gb the row where n lives, i.e. the rank of the party that got bigger
for i from 1 to k while T[i][-1]<>n do od:
gb:=i:

T1:=[op(1..gb-1,T),[op(1..nops(T[gb])-1,T[gb])],op(gb+1..nops(T),T)]:

V:=VoteSeqsC(T1):
{seq(op(ExtVo(k,v,gb)),v in V)}:
end:


#LaneSeqsC(T): a clever (recursive) way of computing LaneSeqs(T).
LaneSeqsC:=proc(T) local k,i,n,T1,gb,V,v,i1:
option remember:
k:=nops(T):

if T=[[]$k] then
 RETURN({[]}):
fi:
n:=add(nops(T[i]),i=1..k):

#gb the row where n lives
for i from 1 to k while T[i][-1]<>n do od:
gb:=i:

T1:=[op(1..gb-1,T),[op(1..nops(T[gb])-1,T[gb])],op(gb+1..nops(T),T)]:

V:=LaneSeqsC(T1):

if gb=k or nops(T[gb])-nops(T[gb+1])>=2 then
RETURN({seq([op(v),gb],v in V)}):
fi:

for i from gb+1 to k while nops(T[i])=nops(T[gb])-1 do od:
i:=i-1:

{seq(seq([op(v),i1],i1=gb..i),v in V)}


end:

RatKc:=proc(T) local n,k,gb,i:
k:=nops(T):
n:=add(nops(T[i]),i=1..k):
for i from 1 to k while T[i][-1]<>n do od:
gb:=i:
for i from gb+1 to k while nops(T[i])=nops(T[gb])-1 do od:
i:=i-gb:
end:






#LaneRange(k,v): all the possible continuations of Lane preferance sequence v, with the last letter chopped (call it v') (with k lanes). Try:
#LaneRange(3,[1,1,2,2,3,3]);
LaneRange:=proc(k,v) local n,T,i,gb,i1:
n:=nops(v):
T:=kLane(k,v):

for i from 1 to k while T[i][-1]<>n do od:
gb:=i:

if gb=k or nops(T[gb])-nops(T[gb+1])>=2 then
RETURN([gb]):
fi:

for i from gb+1 to k while nops(T[i])=nops(T[gb])-1 do od:
i:=i-1:

[seq(i1,i1=gb..i)]:
end:


#VoteRange(k,v): Given a pos. integer k, a word v in {1, ...,k} finds all the ways to extend the chopped v to get the same tableau
VoteRange:=proc(k,v) local T,L,i,i1,n,g:
n:=nops(v):
T:=kVote(k,v):
for i from 1 to k while T[i][-1]<>n do od:
g:=i:

L:=Rank1(SV(v,k)):
 for i from 1 to nops(L) while add(nops(L[i1]),i1=1..i)<g do od:
 convert(L[i],list):

end: