# restart; read `\G:/My Drive/Aek/Card2.1/ProbGen.txt`;

# First version: November 26, 2024

# This program accompanies
# the paper "Prob. Generating Function"
# by Tipaluck Krityakierne
# and myself

with(combinat,choose):
with(LinearAlgebra):


#################################
# Section 1: Data of n cards 
# k shufflings with formulas
#
# Functions:
# MyM(n,k), A1(i,j,n),
# MyP(n,k), Ak(i,j,n,k), 
#
#################################

# Formula of matrix M
# Input: positive integers n and k
# Output: Matrix using
# formulas
# Try: 
# MyM(4,1); 
MyM := proc(n,k) option remember;
local i,j,M;

M := [seq([seq( 2^n*A1(i,j,n),j=1..n)],i=1..n)];
Matrix(M)^k;
end:


# Input: positive integers i,j,n
# Output: formula of the entry 
# (i,j) of matrix of size n
A1 := proc(i,j,n) option remember;
binomial(i-1,j-1)/2^i
+binomial(n-i,j-i)/2^(n-i+1);

end:

# Formula of matrix M
# Input: positive integers n and k
# Output: Matrix using
# formulas
# Try: 
# MyP(4,1); 
MyP := proc(n,k) option remember;
local i,j,P;

P := [seq([seq( A1(i,j,n),j=1..n)],i=1..n)];
Matrix(P)^k;
end:

# Input: positive integers i,j,n,k
# Output:  the entry (i,j) of n cards after
# k shuffles (through the nested sum formula)
# Try:
# Ak(1,1,4,2)
Ak := proc(i,j,n,k) option remember;
local s;

if k=1 then
	return(A1(i,j,n));
fi:
	
add(A1(i,s,n)*Ak(s,j,n,k-1),s=1..n);
end:


#################################
# Section 2: Generating functions
#
# Functions:
# B1(c,i,j,n), MyPL(L,i,j,n),
# EQ5(i,j,n,k), Vec(k),
# 
# GenF(L,i,n,x), 
# FormuF(L,i,n,x), Ind1(L), 
#
#################################

B1 := proc(c,i,j,n) option remember;

if c=1 then
	1/2^(i-1)*binomial(i-1,j-1);
elif c=2 then
	1/2^(n-i)*binomial(n-i,j-i);
else
	ERROR("NoCase");
fi:	

end:

# PL as defined in section2, eq6
# Try: MyPL([1,1],i,j,n);
MyPL := proc(L,i,j,n) option remember;
local s,M;

if nops(L) = 0 then
	if i=j then
		return(1);
	else
		return(0);
	fi:	
fi:
M := [op(2..nops(L),L)]; 
add( B1(L[1],i,s,n)*MyPL(M,s,j,n),s=1..n);

end:


# Value of A_{i,j}^(k) through P_L
# Try:
# EQ5(1,1,4,2);
EQ5 := proc(i,j,n,k) option remember;
local L;

1/2^k*add(MyPL(L,i,j,n), L in Vec(k));

end:

# [0,1] vector of length k
# Try: Vec(2);
Vec := proc(k) option remember;
local L;

if k=0 then
	return({[]});
fi:	

{seq( [1,op(L)], L in Vec(k-1))} union
{seq( [2,op(L)], L in Vec(k-1))};
end:

# Generating function Eq(7)
# Try
# GenF([1,1],2,7,x);
GenF := proc(L,i,n,x) option remember;
local j,F;

F := add( MyPL(L,i,j,n)*x^j,j=1..n);
factor(F);
end:


# Generating function by formula
# Try:
# FormuF([2,1,2],5,7,x);
# GenF([2,1,2],5,7,x);
FormuF := proc(L,i,n,x) option remember;
local g,a,b,t,k;

g := (a,b) -> (a*x+(1-a))^(n-i)
*(b*x+(1-b))^(i-1)*x;

t := Ind1(L);
k := nops(L);

factor(g(t/2^k,(t+1)/2^k));
end:

#Index of L as given in figure 2 and theorem 1
# Try: Ind1([1,1,1]); Ind1([2,1,2]);
Ind1 := proc(L) option remember;
local i,k;

k := nops(L);

add( (L[i]-1)*2^(i-1),i=1..k);
end:


#################################
# Section 3: Applications
#
# Functions:
# Gi(i,n,k,x),
#
#
##################################

# Corollary 3 in the paper
# Try:
# A:=Gi(11,17,7,x):
# B:=add(Ak(11,j,17,7)*x^j,j=1..17);
# simplify(A-B);
Gi := proc(i,n,k,x) option remember;
local j,t;

1/2^k*add( (t/2^k*x+1-t/2^k)^(n-i)
*((t+1)/2^k*x+1-(t+1)/2^k)^(i-1)*x  
,t=0..2^k-1);

end: