# restart; read `\c:/Users/Asus/Google Drive/Aek/
# Card2/NoFeedBack.txt`;

# First version: October 10, 2021
# Second version: January 29, 2022
# Third version: April 18, 2022

# This program accompanies
# the paper "No Feedback? No Worries!
# The art of guessing the right card"
# by Tipaluck Krityakierne
# and myself

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


#################################
# Section 1: Data of n cards 
# k shufflings both empirical 
# and with formulas
#
# Functions:
# Decks(n,k), ShufDeck(D), 
# CountMul(S), PosA(b,S), 
# EmM(n,k),
#
# MyM(n,k), A1(i,j,n),
# BestGuess(n,k), 
#
#################################

##########################
# Section 1.1: Empirical
##########################

# Input: positive integer n and k
# Output: the set of deck of n cards 
# after riffle shuffling k times
# Try: 
# Decks(3,1);
# Decks(4,4);
Decks := proc(n,k) option remember;
local i,j,a,p,A,S,T;

S := [[seq(i,i=1..n)]];
T := [1];

for i from 1 to k do
	A := [0$nops(T)];
for j from 1 to nops(S) do
for a in ShufDeck(S[j]) do
	p := PosA(a[1],S);
	if p > 0 then
		A[p] := A[p]+a[2]*T[j];
	else
		S := [op(S),a[1]];
		A := [op(A),a[2]*T[j]];
	fi:		
od: od: 
	T := A;
od:

{seq( [S[i],T[i]],i=1..nops(S))};	
end:


# Input: Deck D as a list
# Output: the set of all 
# possibilities after one
# riffle shuffling of D.
# Try:
# ShufDeck([1,2,3]);
ShufDeck := proc(D) option remember;
local i,c,t,n,A,S,T;

n := nops(D);
T := [];

for t from 0 to n do
for S in choose(n,t) do
	A := [0$n];
	for i from 1 to nops(S) do
		A[S[i]] := D[i];
	od:	
	c := t+1;
	for i from 1 to n do
		if A[i] = 0 then
			A[i] := D[c];
			c := c+1;
		fi:	
	od:	
	T := [op(T), A];
od: od:
CountMul(sort(T));	

end:


# Input: List of possibilities, S
# Output: Multiplicity of 
# each deck in the list S
# Try: 
# S:=[[1,2],[1,2],[1,2],[2,1]]: 
# CountMul(S);
CountMul := proc(S) option remember;
local i,T,c;

T := {};
c := 1;

for i from 2 to nops(S) do
	if S[i]=S[i-1] then
		c := c+1;
	else
		T := T union {[S[i-1],c]};
		c := 1;
	fi:
od:	

T union {[S[nops(S)],c]};
end:


# Input: positive integer b 
# and the list S
# Output: the position of b 
# in the list S, if no b then
# return 0
# Try: 
# PosA(3,[1,4,2,3,5]);
PosA := proc(b,S) option remember;
local i;

for i from 1 to nops(S) do
	if b = S[i] then
		return(i);
	fi:	
od:
return(0);
end:



# Empirical matrix of size n-by-n
# Input: positive integers n and k
# Output: the list of size n-by-n
# on frequency of number j 
# in position i after k-time 
# riffle shuffles
# Try: 
# EmM(2,2);
# EmM(4,1); 
EmM := proc(n,k) option remember;
local i,s,A;

A := [[0$n]$n];

for s in Decks(n,k) do
for i from 1 to nops(s[1]) do
	A[i][s[1][i]] := A[i][s[1][i]]+s[2];
od: od:

Matrix(A);
end:

##########################
# Section 1.2: Formulas
##########################

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

M := [seq([seq( 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 EmM of size n
A1 := proc(i,j,n) option remember;
2^(n-i)*binomial(i-1,j-1)
+2^(i-1)*binomial(n-i,j-i)

end:


######################################
# Input: positive integer n and k
# Output: list of numbers 
# as the most occuring in each
# position for k-shuffling of n cards.
# Try: 
# BestGuess(3,2);
# [seq(BestGuess(34,2)[i],i=1..34)];
BestGuess := proc(n,k) option remember;
local i,j,m,M,S;

S := Array([{}$n]);
M := MyM(n,k);

for i from 1 to n do
	m := max(M[i]); 
	for j from 1 to n do
		if M[i][j] = m then
			S[i] := S[i] union {j};
		fi:
	od:
od:

S;
end:



###########################
# Section 2: 
# Prove best guess 
# and a formula of E[x_k] 
# 
# Functions:
# A2(i,j,n), P1(i,j,n),
# P2(i,j,n), P3(i,j,n), 
# P4(i,j,n), 
#
###########################

# Value of A2(i,j,n) which 
# is the entry in Matrix P^2
# Try: 
# A2(4,5,9); MyM(9,2)[4,5]/2^(18);
A2 := proc(i,j,n) option remember;

1/4*(P1(i,j,n)+P2(i,j,n)+P3(i,j,n)+P4(i,j,n));
end:

# p[1,1] as in section 2 of the paper
# Try: P1(i,3,n);
P1 := proc(i,j,n) option remember;
local s;

sum(1/2^(i-1)*binomial(i-1,s-1)*
1/2^(s-1)*binomial(s-1,j-1)
,s=1..n);
end:

# p[1,2] as in section 2 of the paper
# Try: 
# n:=12:i:=7: seq(P2(i,j,n),j=1..n);
P2 := proc(i,j,n) option remember;
local s;
sum(1/2^(n-i)*binomial(n-i,s-i)*
1/2^(s-1)*binomial(s-1,j-1)
,s=j..n);

end:

# p[1,3] as in section 2 of the paper
# Try: P3(i,3,n);
P3 := proc(i,j,n) option remember;
local s; 
sum(1/2^(i-1)*binomial(i-1,s-1)*
1/2^(n-s)*binomial(n-s,j-s)
,s=1..j);
end:

# p[1,4] as in section 2 of the paper
# Try: P4(i,3,n);
P4 := proc(i,j,n) option remember;
local s;

sum(1/2^(n-i)*binomial(n-i,s-i)*
1/2^(n-s)*binomial(n-s,j-s)
,s=1..n);

# 4^(i-1)*binomial(n-i,j-i)*3^(j-i);
end:


###################################
# Section 2.1: 
# Split entries for Ak
#
# Functions:
# PL(i,j,n,L), MatPL(n,L),
# A3(i,j,n),
# FL(L,i,n,x),
#
####################################

# Generalize p_L of each a^(k)_{i,j}
# Input: positive integers i,j,n and
# the list L (with entries only 1 or 2)
# Output: p_L according to 
# definition in section 2
# Try:
# PL(5,2,9,[1,1]); P1(5,2,9);
# PL(5,2,9,[1,1,2]); 
# add(PL(5,j,9,[1,1,2]),j=1..9);
PL := proc(i,j,n,L) option remember;
local s,N,T;

if L=[1] then
	return(1/2^(i-1)*binomial(i-1,j-1));
elif L=[2] then
	return(1/2^(n-i)*binomial(n-i,j-i));
fi:

N := nops(L);
T := [op(1..N-1,L)];	

if L[N]=1 then
	add(PL(i,s,n,T)*1/2^(s-1)
	*binomial(s-1,j-1),s=1..n);
elif L[N]=2 then
	add(PL(i,s,n,T)*1/2^(n-s)
	*binomial(n-s,j-s),s=1..n);
else
	ERROR("BadInput");
fi:
end:

# Input: The size of deck n and 
# the list L (of length k) composed
# of only 1 or 2
# Output: Matrix with entries in p_L
# Try:
# MatPL(4,[1,2,2]);
MatPL := proc(n,L) option remember;
local i,j,P;

P := [seq([seq( PL(i,j,n,L),j=1..n)],i=1..n)];
Matrix(P);
end:


# Check the entries of A^(3)(i,j)
# with the actual formula
# Input: positive integers i,j,n
# Output: a^(3)(i,j,n)  
# i.e. entries after 3 shufflings
# Try:
# A3(6,5,6); MyM(6,3)[6,5]/2^(18);
A3 := proc(i,j,n) option remember;
local s,S;

S := {[1,1,1],[1,1,2],[1,2,1],[1,2,2]
,[2,1,1],[2,1,2],[2,2,1],[2,2,2]};

1/8*add(PL(i,j,n,s),s in S);
end:


# Generating function F_L(x) 
# of p_L
# Input: the list L composed of
# entries only 1 or 2, 
# variables i,j,n and x.
# Output: Generating function F_L(x)
# as defined in section 2 of the paper 
# Try: 
# FL([1,2],i,n,x);
# simplify(subs(x=1,FL([1,2,1,2,2],i,n,x)));
FL := proc(L,i,n,x) option remember;
local s,T,F;

if L = [] then
	return(x^i);
fi:	

T := [op(2..nops(L),L)];

if L[1] = 1 then	
	F := sum(1/2^(i-1)*binomial(i-1,s-1)
		*FL(T,s,n,x),s=1..i);
elif L[1] = 2 then
	F := sum(1/2^(n-i)*binomial(n-i,s-i)
		*FL(T,s,n,x),s=i..n);
else
	ERROR();
fi:

return(simplify(factor(F)));	
end:




###################################
# Section 3: 
# Test the approximation of E[X]
# (last theorem)
#
# Functions:
# For1k(k,i,j), Error1k(k,i),
#
###################################

# Test the order of the error 
# in the middle approximation part

# Input: Positive integers k,i,j
# Output: Formula for p_[1,1,...,1](i,j);
# Try: 
# For1k(2,5,3); P1(5,3,10);
For1k := proc(k,i,j) option remember;
binomial(i-1,j-1)*(2^k-1)^(i-j)/2^(k*i)*2^k;
end:


# Calculate the error from normal 
# approximation and the actual
# Input: Positive integers k,i
# Output: the error at the point
# of estimation from normal 
# approximation and the actual
# Try:
# Error1k(3,500);
# i:=16000: evalf(Error1k(4,i)*i^(3/2));
# i:=160000:evalf(Error1k(4,i)*i^(3/2));
# It looks like the error of each term
# is of order 1/i^(3/2)
# seq(evalf(Error1k(5,i)*i^(3/2)),i=32000..32010);
# seq(evalf(Error1k(5,i)*i^(3/2)),i=320000..320010);
Error1k := proc(k,i) option remember;
local j,A,B;

A := 1/sqrt((i-1)*(2^k-1)*2*Pi);
j := floor(i/2^k)+1;
B := For1k(k,i,j)/2^k;
evalf(A-B); 
end: