# restart; read `G:/My Drive/Aek/Card2.2/OptGuess.txt`;

# First version: October 11, 2025
# Second version: September 11, 2026

# My effort to prove Card2 again (first verison)
# Check the pdf from Claude (second version).

#####################################
# Section 1: Introduction from 
# section 1 and 2
#
# Functions: 
# GenG(i,n,k,x), Aij(i,n,k),
# OptGuess(i,n,k), Thm11(i,n,k),
# OptMean(n,k), Thm12(n,k),
#
# FLt(i,n,k,t), Eq23(i,n,k,t),
# 
#####################################

# Main generating function, eq 2.1
# Try:
# A:= GenG(11,17,7,x):

GenG := proc(i,n,k,x) option remember;
local 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:

# Input: position i, number of 
# card n and number of shuffles k
# Output: List of probability a_{i,j}^(k), 1<=j<=n,
# calculate from generating function
# Try: 
# Aij(4,10,3);
# sort(expand(GenG(4,10,3,x)));
# Plot coeff of x^j from G 
# A:=Aij(25,150,3):
# plot([seq(j,j=1..150)],A); 
Aij := proc(i,n,k) option remember;
local j,x,G,A;

G := expand(GenG(i,n,k,x));
A := [seq(coeff(G,x,j),j=1..n)];
return(A);
end:

# Empirical optimal guess from generating function
# Input: position i, number of 
# card n and number of shuffles k
# Output: list of index j and coeff 
# that gives the max. coefficient
# Try: 
# OptGuess(25,150,3);
# seq(OptGuess(i,10,3),i=1..10);
# seq(OptGuess(i,40,3),i=1..40); 
OptGuess := proc(i,n,k) option remember;
local j,A,S,champ;

A := Aij(i,n,k);

S := {};
champ := max(A);

for j from 1 to n do
	if A[j]=champ then
		S := S union {j};
	fi:
od:
[S,champ];	
end:

# From Theorem 1.1: optimal guessing 
# strategy for fixed k
# Input: position i, number of 
# card n and number of shuffles k
# Output: index j that gives 
# the max. coefficient
# Try:
# A:=[seq(op(OptGuess(i,150,3)[1]),i=1..150)];
# B:=[seq(Thm11(i,150,3),i=1..150)];
# evalb(A=B);
Thm11 := proc(i,n,k) option remember;

if i<= floor((n+1)/2) then
	return( floor(i/2^k)+1);
else
	return( n-floor((n+1-i)/2^k));
fi:	

end:


# Empirical optimal guess from generating function
# Input: number of card n and 
# number of shuffles k
# Output: sum of the optimal 
# probability for all position i.
# Try: 
# OptMean(150,3); 
OptMean := proc(n,k) option remember;
local i,A;

A := add( OptGuess(i,n,k)[2],i=1..n);
evalf(A);
end:

# From Theorem 1.2: Expected number of correct 
# guessings from optimal guessing strategy. 
# Input: number of card n and number of shuffles k
# Output: expected value from optimal guessing.
# Try:
# Thm12(150,3);
# OptMean(150,3);
# [seq( OptMean(n,3)-Thm12(n,3), n in {50,100,150,200})]; 
# evalf(-9/(3*7));
Thm12 := proc(n,k) option remember;
local A;

A := 2*sqrt(n)/sqrt((2^k-1)*Pi);
evalf(A);
end:

#################################################
# Input: position i, number of card n, 
# number of shuffles k, component t
# Output: coeff of each j of component t
# Try:
# FLt(25,150,3,3);
# plot the probability distribution of j
# A1:=FLt(25,150,3,1):
# A2:=FLt(25,150,3,2): 
# A3:=FLt(25,150,3,3):
# R1 := plot([seq(j,j=1..150)],A1):
# R2 := plot([seq(j,j=1..150)],A2,color=blue):
# R3 := plot([seq(j,j=1..150)],A3,color=green):
# plots[display]({R1,R2,R3});
FLt := proc(i,n,k,t) option remember;
local x,T,A;

T := x*(t/2^k*x+1-t/2^k)^(n-i)
*((t+1)/2^k*x+1-(t+1)/2^k)^(i-1);

A := [seq(coeff(expand(T),x,j),j=1..n)];
return(A);
end:


# Input: position i, number of card n, 
# number of shuffles k, component t
# Output: prob. of each j in component t
# via prob. distribution of equation 2.3
# Try:
# A:=Eq23(25,150,3,3):
# B:=FLt(25,150,3,3):
# A-B; 
Eq23 := proc(i,n,k,t) option remember;
local s,C1,C2,T,A;

C1 := t/2^k;
C2 := (t+1)/2^k;
T := j-> add( binomial(n-i,s)*C1^s*(1-C1)^(n-i-s)
*binomial(i-1,j-s)*C2^(j-s)*(1-C2)^(i-1-(j-s)),s=0..j);

A := [seq( T(j),j=0..n-1)];
return(A);
end:


######################################
# Section 2: Lemma 3.1,3.2,3.3,3.4 
# from the paper
#
# Functions:
# MyMode(N,k), Lem311(N,k), 
# Lem312(N,k), Lem313(N,k),
# 
# Lem32(i,n,k,t), 
# Lem33(i,j,n,k), CloseT(i,j,n,k),
# Lem34(p,Q),
#
######################################

# Input: position N and number of shuffle k
# Output: value of the Mode for given N and k
# Try:
# MyMode(50,3);
MyMode := proc(N,k) option remember;
local p,x;

p := 1/2^k;
x := floor((N+1)*p);
evalf(binomial(N,x)*p^x*(1-p)^(N-x));
end:

# Test Lemma 3.1(i) on page 4 numerically
# Input: position N and number of shuffle k
# Output: Percentages of the difference and
# percentage of the lower bound.
# Try:
# Lem311(50,3);
Lem311 := proc(N,k) option remember;
local p,x,mode,P,b1,b2;

p := 1/2^k;
mode := {floor((N+1)*p)};

if type((N+1)*p,integer) then
	mode := mode union {(N+1)*p-1};
fi:	

P := MyMode(N,k);

x := min(mode)-1;
b1 := binomial(N,x)*p^x*(1-p)^(N-x);

x := max(mode)+1;
b2 := binomial(N,x)*p^x*(1-p)^(N-x);

evalf([(P-b1)/P,(P-b2)/P, 1/2^k/(N+1)^2/P]);
end:

# Test Lemma 3.1(ii) on page 4 numerically
# Input: position N and number of shuffle k
# Output: Difference between the modes, a lower bound.
# Try:
# Lem312(50,3);
# seq(Lem312(n,3),n=100..120);
Lem312 := proc(N,k) option remember;
local m,p;

m := 2^k;
p := 1/m;

if type((N+1)*p,integer) then
	if MyMode(N,k) = MyMode(N+1,k) then
		return(true);
	else
		return(false,MyMode(N,k),MyMode(N+1,k));
	fi:	
fi:

evalf([MyMode(N,k)-MyMode(N+1,k), 1/m^2/(N+2)^2]);
end:


# Input: position N and number of shuffle k
# Output: Ratio of the leading term 
# and the actual mode.
# Try: 
# Lem313(50,3);
# Lem313(1500000,10);
Lem313 := proc(N,k) option remember;
local p,A;

p := 1/2^k;
A := 1/sqrt(2*Pi*p*(1-p)*N);
[evalf(A)/MyMode(N,k)];
end:


# Input: position i, number of card n, 
# number of shuffles k and curve t
# Output: value of mode H_{i-1} and difference*sqrt(n). 
# Try: Lem32(30,100,3,3);
Lem32 := proc(i,n,k,t) option remember;
local A,B;

if i> floor((n+1)/2) or t<1 or t>2^k-2 then
	ERROR("BadInput");
fi:	

A := max(FLt(i,n,k,t));
B := MyMode(i-1,k);

evalf([B, (B-A)*sqrt(n)]);
end:


# Input: position i, card value j,
# number of card n and number of shuffles k. 
# Output: Values of tails of other curves to j.
# Try:
# Lem33(25,8,150,3);
# Lem33(30,43,150,3);
Lem33 := proc(i,j,n,k) option remember;
local t,m,s,S,eps,actual;

m := 2^k;
s := CloseT(i,j,n,k);
S := {seq(t,t=0..m-1)} minus {s};

actual := add( FLt(i,n,k,t)[j],t in S);
eps := 2*(m-1)*exp(-(n-1)/(2*m^2));
#print(evalf(FLt(i,n,k,s)[j]));

evalf([actual,eps]);
end:



# Input: position i, card value j,
# number of card n and number of shuffles k. 
# Output: Index t of the closest of j to the 
# peak of curve t.
# Try:
# CloseT(30,43,150,3);
CloseT := proc(i,j,n,k) option remember;
local t,s,S,dis;

S := [seq( 1+((n-1)*t+i-1)/2^k,t=0..2^k-1)];
#print(evalf(S));

dis := n;
for t from 0 to 2^k-1 do
	if abs(j-S[t+1]) < dis then
		s := t;
		dis := abs(j-S[t+1]);
	fi:
od:

s;
end:

# Test the claimed statement
# of Lemma 3.4
# Input: probability p and list 
# of probabilities q_l.
# Output: List of differences of 
# the claim from r=1 to D.
# Try:
# Lem34(1/4,[1/2,1/2,2/3]);
Lem34 := proc(p,Q) option remember;
local i,x,D,P,C;

if min(Q) <= p then
	ERROR("BadInput");
fi:	

D := nops(Q);
P := mul( 1-Q[i]+Q[i]*x ,i=1..D);
C := i-> coeff(P,x,i);

print([seq(C(i),i=0..D)]);

[seq( (1-p)*i*C(i)-p*(D-i+1)*C(i-1) ,i=1..D)];
end:


######################################
# Section 3: Cor4.1, Prop5.1 
# and Proofs in section 3 
#
# Functions:
# Cor41(i,n,k), Prop51(i,n,k),
#
#
######################################

# From Corollary 4.1: All optimal guessing 
# strategies for fixed k >=1
# Input: position i, number of 
# card n and number of shuffles k
# Output: set of index j that gives 
# the max. coefficient
# Try:
# A:=[seq(OptGuess(i,151,3)[1],i=1..150)];
# B:=[seq(Cor41(i,151,3),i=1..150)];
# evalb(A=B);
# A:=[seq(OptGuess(i,15,1)[1],i=1..15)];
# B:=[seq(Cor41(i,15,1),i=1..15)];
# evalb(A=B);
# A:=[seq(OptGuess(i,30,1)[1],i=1..30)];
# B:=[seq(Cor41(i,30,1),i=1..30)];
# evalb(A=B);
Cor41 := proc(i,n,k) option remember;

if k >=2 then
	if n mod 2 = 1 and i=(n+1)/2 then
		{floor(i/2^k)+1,n-floor((n+1-i)/2^k)};
	else
		{Thm11(i,n,k)};
	fi:	
#elif k=1 then
#	if n mod 2 = 1 and i=(n+1)/2 then 
#		{floor(i/2^k)+1,floor(i/2^k),n-floor((n+1-i)/2^k)
#		,n-floor((n+1-i)/2^k)+1};	
#	elif i<= floor((n+1)/2) and i>=4 and i mod 2 =0 then
#		{Thm11(i,n,k),Thm11(i,n,k)-1};
#	elif i> floor((n+1)/2) and i<=n-3 and (n-i)mod 2=1 then
#		{Thm11(i,n,k),Thm11(i,n,k)+1};
#	else
#		{Thm11(i,n,k)};
#	fi:	
else
	ERROR("NoCase");
fi:		
end:


# Input: position i, number of card n, 
# number of shuffles k, component t
# Output: bounds and value of a_{i,j}^(k) itself.
# Try:
# Prop51(30,150,3);
# Prop51(30,1500,3);
# Prop51(30,3000,3);
Prop51 := proc(i,n,k) option remember;
local m,eps,A,L,R;

m := 2^k;
eps := 2*(m-1)*exp(-(n-1)/(2*m^2));

A := max(Aij(i,n,k));
L := MyMode(i-1,k)/m;
R := (MyMode(i-1,k)+eps)/m;

evalf([L,A,R]);
end: