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

# First version: May 11, 2021

# This program accompany 
# the paper "The Card Guessing Game:  
# A generating function approach"
# by Tipaluck Krityakierne
# and myself


with(combinat,choose):


##################################
# Section 0: Undetermined Coeff.
# for Central Binomial
# 
# Functions:
# UCCenBi(K,n), Ty1(pow,K,n),
#
##################################


# Undetermined coeff. for
# central binomial from 
# holonomic relation:
# (4*n-2)*f(n-1)-n*f(n)
# Input: number of terms K 
# and symbolic n
# Output: the series approximation
# of binomial(2*n,n)
# Try: 
# UCCenBi(5,n);
# CenBi(5,n);
UCCenBi := proc(K,n) option remember;
local i,a,A,B,T,eq,var;

A := (n^(-1/2)+add(a[i]*n^(-1/2-i),i=1..K))
*4^n/sqrt(Pi);
B := (Ty1(-1/2,K+2,n)
+add(a[i]*Ty1(-1/2-i,K+2-i,n),i=1..K))
*4^(n-1)/sqrt(Pi);

T := simplify(((4*n-2)*B-n*A)/4^n); 
T := numer(T); 
eq := {seq(coeff(T,n,i)=0,i=1..K)};
var := {seq(a[i],i=1..K)};
var := solve(eq,var);
subs(var,A);
end:


# Input: power pow, number 
# of K term, symbolic n
# Output: Taylor series expansion 
# of (n-1)^pow up to K terms
# Try:
# Ty1(1/2,3,n);
Ty1 := proc(pow,K,n) option remember;
local i,x,A,M;

A := taylor( (n+x)^pow,x=0,K); 
M := add( coeff(A,x,i)*x^i,i=0..K-1);
subs(x=-1,M);
end:



#################################
# Section 1: Empirical results
#
# Functions:
# Decks(n), CorrectD(D),
# EmAvg(n), EmVar(n),
#
#################################

# Input: positive integer n
# Output: the list of
# possible decks of n cards 
# after 1-time riffle shuffle
# Try: 
# Decks(3);
# Decks(4);
Decks := proc(n) option remember;
local i,t,T,S,D1,A,c;

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]] := i;
	od:	
	D1 := {seq(i,i=1..n)} minus {op(S)};
	D1 := sort([op(D1)]); 
	c := t+1;
	for i from 1 to nops(D1) do
		A[D1[i]] := c;
		c := c+1;
	od:	
	T := [op(T), A]; 
od: od:
sort(T);
end:



# Perform a guessing on the deck
# of cards.
# Input: Deck D
# Output: Number of correct guesses
# according to the optimal strategy
# Try:
# CorrectD([1,2,3,4]);
# CorrectD([2,3,1,4]);
CorrectD := proc(D) option remember;
local i,n,correct,m,k,C;

i := 1:
n := nops(D);

while i <= n and D[i] = i do
	i := i+1;
od:

correct := i-1;

if i = n+1 then
	return(correct);
fi:	

m := i;
k := D[i];
C := [k-m,n-k];
i := i+1;

while i <= n do
	if (C[1] >= C[2] and D[i] < k) or 
	   (C[1] < C[2] and D[i] > k)	then
		correct := correct+1;
	fi:		
	
	if D[i] < k then
		C[1] := C[1] - 1;
	elif D[i] > k then	
		C[2] := C[2] - 1;	
	else	
		ERROR();
	fi:		
	i := i+1;
od:
#print(C);
correct;
end:


# Input: positive integer n
# Output: average number of 
# correct guesses of the deck 
# of n cards computed empirically
# Try: EmAvg(3);
EmAvg := proc(n) option remember;
local s,S;

S := Decks(n);
add( CorrectD(s),s in S)/nops(S);
end:


# Input: positive integer n
# Output: Variance of the number of 
# correct guesses of the deck 
# of n cards computed empirically
# Try: EmVar(3);
EmVar := proc(n) option remember;
local s,S;

S:=Decks(n);
add( (CorrectD(s)-EmAvg(n))^2,s in S)/nops(S);
end:


#################################
# Section 2: Generating function
# method for higher-order moments
#
# Functions:
# EmGenD(n,q), GenD(n,q),
# RecGenG(n,q), RecGenF(a,b,q), 
# ForGenF(a,b,q), ForGenG(n,q), 
#
#################################


# Empirical generating function 
# for number of correct guesses 
# in the deck of n cards
# Input: positive integer n 
# and symbolic q
# Try: EmGenD(2,q);
EmGenD := proc(n,q) option remember;
local s,D;

D := add( q^CorrectD(s),s in Decks(n));
expand(D);
end:


# Generating function of Deck
# of n cards calculated by
# recurrences
# Input: positive integer n 
# and symbolic q
# Try: 
# [seq(GenD(n,q)-EmGenD(n,q),n=1..12)];
GenD := proc(n,q) option remember;
local i,D;

if n=0 then
	return(1);
fi:

D := q*GenD(n-1,q)+ RecGenG(n,q);
expand(D);
end:


# Generating function G
# as mentioned in the paper
# Input: positive integer n 
# and symbolic q
# Try: RecGenG(5,q);
RecGenG := proc(n,q) option remember;
local i,G;
G := q^n+add(RecGenF(n-1-i,i,q),i=1..n-1);
expand(G);
end:


# Generating function F 
# from the simple recurrence
# Input: non-negative integers 
# a,b and symbolic q
# Try: 
# seq(RecGenF(n,n,q),n=1..8);
RecGenF := proc(a,b,q) option remember;
local m,n,F;
m := max(a,b);
n := min(a,b);

if n=0 then
	q^m;
else
	F := q*RecGenF(m-1,n,q)+RecGenF(m,n-1,q);
	expand(F);
fi:	
end:

# Generating function F 
# from the formula from 
# proposition in the paper
# Input: non-negative integers 
# a,b and symbolic q
# Try: 
# [seq(ForGenF(n,1,q),n=1..10)];
# [seq(ForGenF(n,n-1,q)-RecGenF(n,n-1,q),n=1..10)];
ForGenF := proc(a,b,q) option remember;
local i,m,n,A;

m := max(a,b);
n := min(a,b);

A := add( q^(m+n-i)*(binomial(m+n,i)
-binomial(m+n,i-1)) ,i=0..n);

expand(A);
end:

# Generating function G
# from the formula from 
# section 2.3.1
# Input: positive integers n
# and symbolic q
# Try: [seq(ForGenG(n,q)-RecGenG(n,q),n=1..10)];
ForGenG := proc(n,q) option remember;
local j,G,k;

if n = 0 then
	return(1);
elif n mod 2 = 0 then
	k := n/2;
	G := q^(2*k)-q^(2*k-1)
	+2*add(ForGenF(2*k-1-j,j,q),j=0..k-1);
else
	k := (n-1)/2;
	G := q^(2*k+1)-q^(2*k)+ForGenF(k,k,q)
	+2*add(ForGenF(2*k-j,j,q),j=0..k-1);
fi:

expand(G);
end :


##################################
# Section 2.2: Evaluate 
# G^(r)_n(q) at q=1.
#
# Functions:
# RecGr1(n,r), ForGr1(n,r),
# Fall(n,r), CloseGr(r,k),
# 
##################################

# Values for G^(r) in Proposition 3
# calculate empirically from RecGenG
# Input: positive integers n
# and non-negative integer r.
# Try: [seq(RecGr1(n,1),n=1..30)];
# r:=3:[seq(ForGr1(n,r)-RecGr1(n,r),n=1..30)];
RecGr1 := proc(n,r) option remember;
local q;

if r = 0 then
	subs(q=1,RecGenG(n,q));
else
	subs(q=1,diff(RecGenG(n,q),q$r));
fi:	
end:


# Formula (8) and (9) 
# in Proposition 3
# Input: positive integers n
# and non-negative integer r.
# Try: [seq(ForGr1(n,4),n=1..30)];
ForGr1 := proc(n,r) option remember;
local i,k,G;

if n mod 2 = 0 then
	k := n/2;
	G := Fall(2*k,r)-Fall(2*k-1,r)
	+2*add( (k-i)*Fall(2*k-1-i,r)*
	(binomial(2*k-1,i)-binomial(2*k-1,i-1))
	,i=0..k-1);
else
	k := (n-1)/2;
	G := Fall(2*k+1,r)-Fall(2*k,r)
	+2*add( (k-i+1/2)*Fall(2*k-i,r)*
	(binomial(2*k,i)-binomial(2*k,i-1))
	,i=0..k);
fi:	

expand(G);
end:


# Input: integer n and 
# non-negative integer r
# Output: Falling factorial (n)_r
# Try: Fall(5,4);
# seq(Fall(n,1),n=0..10);
Fall := proc(n,r) option remember;
local i;
mul(n-i,i=0..r-1);
end:


# Closed form in k for G^(r)_n(q)|{q=1}
# Odd and even cases
# Input: non-negative integer r
# and symbolic k
# Try: CloseGr(1,k);
CloseGr := proc(r,k) option remember;
local i,A,B;

A := Fall(2*k+1,r)-Fall(2*k,r)
	+2*sum( (k-i)*Fall(2*k-i,r)*
	(binomial(2*k,i)-binomial(2*k,i-1))
	,i=0..k)+sum( Fall(2*k-i,r)*
	(binomial(2*k,i)-binomial(2*k,i-1))
	,i=0..k);
B := Fall(2*k,r)-Fall(2*k-1,r)
	+2*sum( (k-i)*Fall(2*k-1-i,r)*
	(binomial(2*k-1,i)-binomial(2*k-1,i-1))
	,i=0..k-1);

[expand(A),expand(B)];
end:




###############################
# Section 2.3: Solve the 
# difference equation by 
# the method of undetermined 
# coeff.
#
# Functions:
# LibGr(r,s,k), CheckGr(n,r),
# CenBiShift(t,k,n), SolveDr(r,s,K,n),
# SolveDrS1(r,s,n), SolveDrS2(r,s,K,n),
# ForDr(r,s,K,n), CheckDr(n,r,K),
# Part2Shift(r,s,t,K,n), Ty2(pow,K,n), 
#
###############################

# Library to keep the formulas
# for G^(r)_{2k+s}(q)|{q=1},
# where s = 0 or 1. 
# (Formula from CloseGr)
# Input: positive integer r,
# s=0 or 1, symbolic k 
# Output: [P1,P2,P3] where
# G = P1*4^k+P2*binomial(2*k,k)+P3
# Try:
# LibGr(1,1,k);
LibGr := proc(r,s,k) option remember;
local T;

if r = 1 then
	if s=1 then 
		T := [(2*k-1)/2,(4*k+1)/2,1];
	elif s=0 then
		T := [(k-1)/2,k,1];
	fi:	
elif r = 2 then
	if s=1 then 
		T := [k^2-k/2+1,4*k^2-3*k-1,4*k];
	elif s=0 then
		T := [(2*k^2-3*k+3)/4,(2*k-3)*k,4*k-2];
	fi:	
elif r = 3 then
	if s=1 then 
		T := [k^3-k/4-3,6*k^3-25*k^2/2+8*k+3
			,12*k^2-6*k];
	elif s=0 then
		T := [k^3/2-3*k^2/4+k/4-3/2
			,3*k^3-10*k^2+10*k,12*k^2-18*k+6];
	fi:		
elif r = 4 then
	if s=1 then 
		T := [k^4+k^3-37/4*k^2+17/4*k+12
		,8*k^4-30*k^3+46*k^2-30*k-12
		,32*k^3-48*k^2+16*k];
	elif s=0 then
		T := [k^4/2-k^3/2-37/8*k^2+55/8*k+15/4
		,4*k^4-22*k^3+46*k^2-40*k
		,32*k^3-96*k^2+88*k-24];
	fi:	
elif r = 5 then
	if s=1 then 
		T := [k^5+5/2*k^4-135/4*k^3+595/8*k^2-233/8*k-60
		,10*k^5-115/2*k^4+142*k^3-417/2*k^2+144*k+60
		,80*k^4-240*k^3+220*k^2-60*k];
	elif s=0 then
		T := [k^5/2-145/8*k^3+255/4*k^2-519/8*k-45/4
		,5*k^5-40*k^4+131*k^3-228*k^2+192*k
		,80*k^4-400*k^3+700*k^2-500*k+120];
	fi:	
elif r = 6 then
	if s=1 then 
		T := [k^6+9/2*k^5-335/4*k^4+2655/8*k^3
		-4333/8*k^2+198*k+360,12*k^6-97*k^5+337*k^4
		-737*k^3+1145*k^2-840*k-360,192*k^5-960*k^4
		+1680*k^3-1200*k^2+288*k];
	elif s=0 then
		T := [k^6/2+3/4*k^5-365/8*k^4+4065/16*k^3
		-4679/8*k^2+8253/16*k+315/8,6*k^6-65*k^5
		+296*k^4-769*k^3+1264*k^2-1092*k,192*k^5
		-1440*k^4+4080*k^3-5400*k^2+3288*k-720];
	fi:		
elif r = 7 then
	if s=1 then 
		T := [k^7+7*k^6-343/2*k^5+1015*k^4-46571/16*k^3
		+66493/16*k^2-11733/8*k-2520,
		14*k^7-301/2*k^6+686*k^5-3859/2*k^4+4325*k^3
		-7445*k^2+5760*k+2520,
		448*k^6-3360*k^5+9520*k^4-12600*k^3
		+7672*k^2-1680*k];
	elif s=0 then
		T := [k^7/2+7/4*k^6-749/8*k^5+11725/16*k^4
		-21539/8*k^3+82117/16*k^2-33471/8*k-315/2
		,7*k^7-98*k^6+581*k^5-1992*k^4+4686*k^3
		-7912*k^2+7248*k,448*k^6-4704*k^5+19600*k^4
		-41160*k^3+45472*k^2-24696*k+5040];
	fi:		
elif r = 8 then
	if s=1 then 
		T := [k^8+10*k^7-623/2*k^6+2527*k^5
		-170191/16*k^4+208985/8*k^3-556901/16*k^2
		+96429/8*k+20160,
		16*k^8-220*k^7+1260*k^6-4208*k^5+11140*k^4
		-28812*k^3+56104*k^2-45360*k-20160
		,1024*k^7-10752*k^6+44800*k^5
		-94080*k^4+103936*k^3-56448*k^2+11520*k];
	elif s=0 then
		T := [k^8/2+3*k^7-679/4*k^6+3507/2*k^5
		-290591/32*k^4+436107/16*k^3-1494453/32*k^2
		+580959/16*k+2835/4
		,8*k^8-140*k^7+1036*k^6-4372*k^5
		+12684*k^4-30000*k^3+55920*k^2
		-55296*k,1024*k^7-14336*k^6+82432*k^5
		-250880*k^4+433216*k^3-420224*k^2+209088*k-40320];
	fi:		
elif r = 9 then
	if s=1 then 
		T := [k^9+27/2*k^8-1041/2*k^7+21987/4*k^6
		-500199/16*k^5+3544443/32*k^4-1999061/8*k^3
		+10261449/32*k^2-1758717/16*k-181440
		,18*k^9-615/2*k^8+2148*k^7-8079*k^6+20982*k^5
		-125799/2*k^4+216828*k^3-481170*k^2+403200*k
		+181440,
		2304*k^8-32256*k^7+185472*k^6
		-564480*k^5+974736*k^4-945504*k^3
		+470448*k^2-90720*k];
	elif s=0 then
		T := [k^9/2+9/2*k^8-1131/4*k^7
		+14805/4*k^6-809151/32*k^5+3389841/32*k^4
		-9031507/32*k^3+14499495/32*k^2
		-5473935/16*k-14175/4
		,9*k^9-192*k^8+1722*k^7-8544*k^6+27393*k^5
		-72816*k^4+200268*k^3-444480*k^2+478080*k
		,2304*k^8-41472*k^7+314496*k^6-1306368*k^5
		+3232656*k^4-4844448*k^3+4252464*k^2
		-1972512*k+362880];
	fi:	
elif r = 10 then
	if s=1 then 
		T := [k^10+35/2*k^9-1635/2*k^8+43395/4*k^7
		-1269807/16*k^6+11979555/32*k^5
		-19223395/16*k^4+82622365/32*k^3
		-3245499*k^2+8850915/8*k+1814400
		,20*k^10-415*k^9+3459*k^8-14082*k^7
		+25878*k^6-29955*k^5+300991*k^4
		-1828268*k^3+4626852*k^2-3991680*k
		-1814400,5120*k^9-92160*k^8+698880*k^7
		-2903040*k^6+7183680*k^5-10765440*k^4
		+9449920*k^3-4383360*k^2+806400*k];
	elif s=0 then
		T := [k^10/2+25/4*k^9-885/2*k^8+57045/8*k^7
		-1971627/32*k^6+21608265/64*k^5
		-79554605/64*k^4+196619765/64*k^3
		-303054453/64*k^2+112354605/32*k+155925/8
		,10*k^10-255*k^9+2712*k^8-15306*k^7
		+48414*k^6-95835*k^5+277808*k^4
		-1355964*k^3+3949296*k^2-4625280*k
		,5120*k^9-115200*k^8+1113600*k^7
		-6048000*k^6+20247360*k^5-43092000*k^4
		+57894400*k^3-46908000*k^2+20531520*k
		-3628800];
	fi:	
elif r = 11 then
	if s=1 then 
		T := [k^11+22*k^10-4895/4*k^9+39765/2*k^8
	-2895057/16*k^7+8672433/8*k^6-292721605/64*k^5 
	+442703635/32*k^4-1849977811/64*k^3 
	+1149514839/32*k^2-48837375/4*k-19958400
	,22*k^11-1089/2*k^10+5324*k^9-22704*k^8+836295/2*k^6
	-1535794*k^5-263271*k^4+17111968*k^3
	-49279548*k^2+43545600*k+19958400
	,11264*k^10-253440*k^9+2449920*k^8-13305600*k^7
	+44544192*k^6-94802400*k^5+127367680*k^4
	-103197600*k^3+45169344*k^2-7983360*k];
	elif s=0 then
		T := [1/2*k^11+33/4*k^10-660*k^9+102465/8*k^8
	-4348707/32*k^7+59828769/64*k^6-285595475/64*k^5
	+966963855/64*k^4-2276853843/64*k^3+427576941/8*k^2
	-627055965/16*k-467775/4
	,11*k^11-330*k^10+4092*k^9-25608*k^8
	+67683*k^7+107706*k^6-1193762*k^5+1760568*k^4
		+8663256*k^3-38949696*k^2+49524480*k
	,11264*k^10-309760*k^9+3717120*k^8
	-25555200*k^7+111072192*k^6-317523360*k^5+601379680*k^4
	-740036000*k^3+561157344*k^2-233830080*k+39916800];
	fi:
elif r = 12 then
	if s=1 then 
		T := [k^12+27*k^11-7051/4*k^10+137445/4*k^9
	-6073617/16*k^8+44781561/16*k^7-949109953/64*k^6
	+3735680355/64*k^5-10901314721/64*k^4+22410433353/64*k^3
	-13838239953/32*k^2+1173982275/8*k+239500800
	,24*k^12-698*k^11+7898*k^10-34254*k^9
	-118800*k^8+2346090*k^7-13527174*k^6
	+38413574*k^5-27384524*k^4-176139432*k^3
	+575605296*k^2-518918400*k-239500800
	,24576*k^11-675840*k^10+8110080*k^9
	-55756800*k^8+242339328*k^7-692778240*k^6
	+1312101120*k^5-1614624000*k^4+1224343296*k^3
	-510174720*k^2+87091200*k];
	elif s=0 then
		T := [1/2*k^12+21/2*k^11-7579/8*k^10+174075/8*k^9
	-8873337/32*k^8+74455623/32*k^7-1770708577/128*k^6
	+7711702845/128*k^5-24698677091/128*k^4+56362385739/128*k^3
	-5208269499/8*k^2+15173342145/32*k+6081075/8
	,12*k^12-418*k^11+5962*k^10-40524*k^9
	+60192*k^8+1190766*k^7-10573566*k^6+41823424*k^5
	-74149624*k^4-39690048*k^3+423241344*k^2-581368320*k
	,24576*k^11-811008*k^10+11827200*k^9-100362240*k^8
	+549001728*k^7-2025644544*k^6+5122381440*k^5-8831180160*k^4
	+10104775296*k^3-7244062848*k^2+2893052160*k-479001600];
	fi:				
else
	ERROR(NoCaseYet);
fi:

T;
end:


# Numerical values of 
# G^(r)_{n}(q)|{q=1},
# from the formula of LibGr
# Input: positive integer n and r
# 
# Try:
# [seq(CheckGr(n,3),n=1..30)];
# [seq(RecGr1(n,3)-CheckGr(n,3),n=1..30)];
# r:=12:[seq(RecGr1(n,r)-CheckGr(n,r),n=1..60)];
CheckGr := proc(n,r) option remember;
local s,k,T;

s := n mod 2;
T := LibGr(r,s,k); 
T := T[1]*4^k+T[2]*binomial(2*k,k)+T[3];
eval(subs(k=(n-s)/2,T));

end:


# Input: shift t, Number of term K, 
# and symbolic n
# Output: Taylor series approximation 
# of binomial(n-t,(n-t)/2) 
# up to K terms
# Remark: has only t=0,1,2 case.
# Try: CenBiShift(1,5,n);
# simplify(CenBiShift(0,5,2*n)/UCCenBi(5,n));
CenBiShift := proc(t,K,n) option remember;
local T;

T := simplify(CenBi(K,(n-t)/2)); 
if t=1 then  
	T := subs((n-1)^(-(2*K+1)/2)
	=Ty1(-(2*K+1)/2,K+1,n),T); 
elif t=2 then  
	T := subs((n-2)^(-(2*K+1)/2)
	=Ty2(-(2*K+1)/2,K+1,n),T);
elif t <> 0 then
	ERROR(NoCaseYet);
fi:
	
T := simplify(expand(T)); 
end:


# Input: Positive integer r,
# case s=1 for odd n and 
# case s=0 for even n,
# the number of term of accuracy K,
# symbolic n
# Output: [P1,P2] where
# P1*4^k+P2*binomial(2*k,k) is
# an approximation of D^(r)_n(q)|{q=1} 
# up to K terms
# Try: SolveDr(10,1,15,n);
SolveDr := proc(r,s,K,n) option remember;

if K < r then
	ERROR("Need K >= r");
fi:

if r=0 then
	return([2^n,0]);
fi:
	
[SolveDrS1(r,s,n),SolveDrS2(r,s,K,n)];
end:


# Solve for the part of 4^k
# Try: SolveDrS1(1,1,n);
# SolveDrS1(2,1,n);
SolveDrS1 := proc(r,s,n) option remember;
local i,k,a,T,A,B,S,eq,var;

if s = 1 then
	T := subs(k=(n-1)/2,(LibGr(r,1,k)[1]
	+LibGr(r,0,k)[1])*4^k); 
elif s = 0 then
	T := subs(k=n/2, LibGr(r,0,k)[1]*4^k
	+LibGr(r,1,k-1)[1]*4^(k-1)); 
else
	ERROR(NoCase);
fi:

T:= simplify(expand(T)
+r*subs(n=n-1,SolveDr(r-1,(s+1) mod 2,r-1,n)[1])
+r*subs(n=n-2,SolveDr(r-1,s,r-1,n)[1])); 
A := add(a[i]*n^i,i=0..r)*2^n; 
B := add(a[i]*(n-2)^i,i=0..r)*2^(n-2);

S := simplify((A-B-T)/2^n); 
S := numer(S);  

eq := {seq(coeff(S,n,i)=0,i=0..r)}; 
var := {seq(a[i],i=0..r)};
var := solve(eq,var); 
factor(subs(var,A));
end:



# Solve for the part of 
# binomial(2*k,k)
# Try: SolveDrS2(1,1,5,n);
SolveDrS2 := proc(r,s,K,n) option remember;
local i,k,a,T,A,B,S,deg,eq,var;

if s = 1 then
	T := subs(k=(n-1)/2,LibGr(r,1,k)[2]
	+LibGr(r,0,k)[2])*CenBiShift(s,K,n);
elif s = 0 then
	T := subs(k=n/2,LibGr(r,0,k)[2]*CenBiShift(0,K,n)
	+LibGr(r,1,k-1)[2]*CenBiShift(2,K,n));
else
	ERROR(NoCase);
fi:

T := simplify(expand(T)
+r*Part2Shift(r-1,(s+1) mod 2,1,K,n)
+r*Part2Shift(r-1,s,2,K,n) ); 

A := add(a[i]*n^(r-1/2-i),i=0..K)*2^n; 
B := add(a[i]*Ty2(r-1/2-i,K+1-i,n),i=0..K)*2^(n-2);

S := simplify((A-B-T)/2^n); 
S := numer(S); 

deg := degree(S,n);
eq := {seq(coeff(S,n,i)=0,i=deg-K..deg)}; 
var := {seq(a[i],i=0..K)};
var := solve(eq,var);
subs(var,A);
end:


# Formula for Dr with K+1 terms
# of accuracy (obtained from
# SolveDr)
# Input: Positive integer r,
# case s=1 for odd n and 
# case s=0 for even n,
# the number of term of accuracy K,
# symbolic n
# Try: ForDr(1,1,5,n);
ForDr := proc(r,s,K,n) option remember;
SolveDr(r,s,K,n)[1]+SolveDr(r,s,K,n)[2];
end:


# Check the numerical result 
# of ForDr
# Input: Positive integers n, r and 
# the number of term of accuracy K
# Try: CheckDr(51,1,5);
# CheckDr(150,1,5);
# CheckDr(150,1,15);
# seq(CheckDr(200,r,15),r=1..10);
CheckDr := proc(n,r,K) option remember;
local M,A,B;

Digits:=50:
A := evalf(subs(M=n,ForDr(r,n mod 2,K,M)));
B := diff(GenD(n,q),q$r);
B := subs(q=1,B);

(A-B)/B;
end:



# Shift for the second part of 
# solution, i.e. binomial(2*k,k)
# Input: Positive integer r,
# case s=1 for odd n and 
# case s=0 for even n, shift t,
# the number of term of accuracy K,
# symbolic n
# Try: Part2Shift(1,1,1,5,n);
Part2Shift := proc(r,s,t,K,n) option remember;
local T,deg,d;

T := simplify(SolveDr(r,s,K,n)[2]); 
deg := degree(simplify(numer(T)/2^n),n);
T := subs(n=n-t,T);
d := -(deg-1/2-(r-1)); 

if t=1 then  
	T := subs((n-1)^d=Ty1(d,K+1,n),T); 
elif t=2 then  
	T := subs((n-2)^d=Ty2(d,K+1,n),T);
elif t <> 0 then
	ERROR(NoCaseYet);
fi:

simplify(expand(T)); 
end:



# Input: power pow, number 
# of K term, symbolic n
# Output: Taylor series expansion 
# of (n-2)^pow up to K terms
# Try:
# Ty2(1/2,3,n);
Ty2 := proc(pow,K,n) option remember;
local i,x,A,M;

A := taylor( (n+x)^pow,x=0,K); 
M := add( coeff(A,x,i)*x^i,i=0..K-1);
subs(x=-2,M);
end:



###############################
# Section 2.4: Formula and 
# check moments about the mean
#
# Functions:
# EmStMo(n,r), EmMoMean(n,r),
# ForFacM(r,s,K,n),
# Stir2(n,k), StMo(r,s,K,n),
# MoMean(r,s,K,n), 
# TestMoMean(n,r,K), 
#
#
################################

# Input: Positive integer n and 
# non-negative integer r 
# Output: Empirical calculation
# of straight moment
# Try: EmStMo(50,1);
EmStMo := proc(n,r) option remember;
local i,q,T;

T := GenD(n,q);
for i from 1 to r do
	T := q*diff(T,q); 
od:

subs(q=1,T)/2^n;
end:


# Input: Positive integer n and 
# non-negative integer r 
# Output: Empirical calculation
# of moment about the mean 
# Try: EmMoMean(50,2);
EmMoMean := proc(n,r) option remember;
local i;

add( (-1)^(r-i)*binomial(r,i)
*EmStMo(n,i)*EmStMo(n,1)^(r-i),i=0..r);
end:


# Input: Positive integer r,
# case s=1 for odd n and 
# case s=0 for even n, 
# the number of term of accuracy K,
# symbolic n
# Outuput: Formula in n of 
# the r^{th} factorial moment up to
# K terms
# Try: ForFacM(1,1,5,n);
ForFacM := proc(r,s,K,n) option remember;
expand(simplify(ForDr(r,s,K,n)/2^n));
end:


# Input: non-negative integers n, k
# Output: Stirling number of second kind
# Try: n:=10:
# factor(add(Stir2(n,k)*Fall(x,k),k=0..n));
Stir2 := proc(n,k) option remember;
local i;

1/k!*add( (-1)^i*binomial(k,i)*(k-i)^n
,i=0..k);
end:


# Input: Positive integer r,
# case s=1 for odd n and 
# case s=0 for even n, 
# the number of term of accuracy K,
# symbolic n
# Output: Formula in n of 
# the r^{th} straight moment up to K 
# (negative power) terms
# Try: StMo(1,1,5,n);
StMo := proc(r,s,K,n) option remember;
local i,A;

A := add(Stir2(r,i)*ForFacM(i,s,K+i,n),i=0..r);
sort(expand(simplify(A)));
end:


# Input: Positive integer r,
# case s=1 for odd n and 
# case s=0 for even n, 
# the number of term of accuracy K,
# symbolic n
# Output: Formula in n of 
# the r^{th} moment about the mean
# up to K terms
# Try: MoMean(2,1,5,n);
MoMean := proc(r,s,K,n) option remember;
local i,A,T;

A := add( (-1)^(r-i)*binomial(r,i)
*StMo(i,s,K+i,n)*StMo(1,s,K+1+r-i,n)^(r-i),i=0..r);
A:= sort(expand(simplify(A)));

T := add(coeff(A,n^(i+1/2))*n^(i+1/2)
,i=-K-1..floor(r/2));
A := A-add(coeff(A,n^(i+1/2))*n^(i+1/2)
,i=-100*K..floor(r/2));
A := expand(simplify(expand(A)));

A:=T + add(coeff(A,n,i)*n^i,i=-K..floor(r/2));
sort(A);
end:


# Test numerically the relative
# difference of the moment about 
# the mean,  formula vs empirical
# Input: Positive integers n,r and K.
# Output: relative difference under
# the r^{th} moment about the mean 
# (correct up to K terms) of X_n
# Try: 
# TestMoMean(150,2,5);
# TestMoMean(150,3,5);
TestMoMean := proc(n,r,K) option remember;
local M,A,B;

Digits := 50:
A := MoMean(r,n mod 2,K,M);
A := subs(M=n,A);
B := EmMoMean(n,r);

evalf((A-B)/B);
end:



################################
# Section 3: k-shuffle
#
# Section 3.1: 
# Empirical k-shuffle
# 
# Functions:
# ShufDeck(D), CountMul(S),
# KShuff(n,k), PosA(b,S),
#
################################

# Input: Deck D as a list
# Output: the list 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 n and k
# Output: Deck of n cards after 
# shuffling k times
# Try: 
# KShuff(3,1);
# KShuff(4,4);
KShuff := 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: positive integer b 
# and the list S
# Output: the position of b 
# in the list S, if no b then
# return 0
# Try: 
# PosA(1,[2,1,3]);
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:




#############################
# Section 3.2: 
# Faster version of section 3.1,
# Creat increasing seq.
# use binomial(n+i,n) for 
# multiplicities, then 
# investigate the best guess
#
# Functions:
# InSeq(n,r), Adjoin(A,inc,n,r),
# Worpitzky(n,k),
# SlowBG(n,k,S), Inc(R), 
#
############################

# Input: positive integer n
# Output: List of set of upto r
# increasing subsequence of 
# permutation of length n
# Try:
# InSeq(4,2);
InSeq := proc(n,r) option remember;
local i;
Adjoin([],0,n,r);
end:


# Input: list of top elements 
# A of whole permutation,
# inc = number of increasing seq in A,
# length of permutation n and k shuffling
# Output: List of sets of permutations 
# with upto r increasing sequences
# that started with A
# Try: 
# Adjoin([],0,5,5);
# Adjoin([4],1,5,4);
Adjoin := proc(A,inc,n,r) option remember;
local i,s,S,T,B,MM;

MM := [{}$r];

if inc > r then
	return(MM);
elif nops(A)=n then
	MM[inc] := {A};
	return(MM);
fi:

S := {seq(i,i=1..n)} minus {op(A)};

for s in S do
	T := [op(A),s];
	if member(s-1,{op(A)}) then
		B := Adjoin(T,inc,n,r);
	elif inc < r then
		B := Adjoin(T,inc+1,n,r);
	else
		B := [{}$r];
	fi:	
	MM := [seq( B[i] union MM[i],i=1..r)];
od:

MM;
end:


# Double check with Worpitzky identities
# Input: positive integer n
# Output: sum from Worpitzky identity
# (suppose to be 2^(k*n))
# Try: seq(ifactor(Worpitzky(n,1)),n=1..8);
Worpitzky := proc(n,k) option remember;
local m,A;

A := InSeq(n,2^k);
add(nops(A[m])*binomial(n+2^k-m,n) ,m=1..2^k);
end:


# Best Guess of the next element 
# Input: n-permuation k shuffles 
# starting with the list S, 
# Output: the next most 
# likely elements after S
# Try: 
# SlowBG(5,2,[3]); 
# SlowBG(10,2,[5]);
SlowBG := proc(n,k,S) option remember;
local i,j,m,R,T,MM;

T := [0$n];

for i from 1 to n do
if member(i,{op(S)}) = false then	
	R := [op(S),i];	
	MM := Adjoin(R,Inc(R),n,2^k);
	T[i] := add(binomial(n+2^k-m,2^k-m)*
			nops(MM[m]),m=1..2^k);
fi:			
od:	

for i from 1 to n do
	if T[i] = max(op(T)) then
		return([i,T]);
	fi:
od:
ERROR("NotExpect");	
end:


# Input: List R
# Output: number of 
# increasing seq. of R
# Try: Inc([3,2,1,4]);
Inc := proc(R) option remember;
local i,j,count;	

count := 0;
for i from 1 to nops(R) do
	if member(R[i]-1,{seq(R[j],j=1..i-1)}) 
	= false then
		count := count+1;
	fi:	
od:

count;		
end:


##################################
# Section 3.3: Faster BestGuess
#
# Function:
# RecInc(n,S), Push(S),
# NumAfter(n,r,M),
# FastBestGuess(n,r,S),
# FastBG(n,k,S),
#
##################################

# Input: positive integer n and 
# the list S
# Output: The distribution of 
# number of increasing seq. (not include S)
# starting with S on perm. of n.
# Try: RecInc(10,[2,5,3,9]);
RecInc := proc(n,S) option remember;
local t,q,r,T,Q,R;

if nops(S) = n then
	if member(1,S) then
		return([1,0$n]);
	else
		return([0,1,0$(n-1)]);
	fi:	
fi:	

T := {op(S)};
Q := {seq(t+1,t in T)} minus {n+1} minus T;
R := {seq(t,t=1..n)} minus Q minus T;

add(RecInc(n,T union {q}), q in Q)+
add(Push(RecInc(n,T union {r})), r in R);
end:


# Input: List S
# Output: Shift the list S by 1
Push := proc(S) option remember;
[0,op(1..nops(S)-1,S)];
end:


# Number of n-card decks that 
# started with S with at most r 
# increasing subseq 
# (for multiplicities)
# Try: 
# [seq(NumAfter(5,4,[3,i]),i=[1,2,4,5])];
# SlowBG(5,2,[3]);
NumAfter := proc(n,r,M) option remember;
local t;
add(binomial(n+r-Inc(M)-t,n)
*RecInc(n,M)[t+1], t=0..min(r-Inc(M),n));
end:


# Input: positive integer n, 
# r is the maximum number 
# of inc. subseq. For example
# 2 shuffles -> r=4, 
# 3 shuffles -> r=8, ...,
# starting with list of permutation S
# Output: the next most 
# likely elements after S
# Try: 
# FastBestGuess(5,4,[3]);
# FastBestGuess(10,4,[5]); 
FastBestGuess := proc(n,r,S) option remember;
local i,T,mx;

T := [0$n];

for i from 1 to n do
if member(i,{op(S)}) = false then
	T[i] := NumAfter(n,r,[op(S),i]) ;
fi:			
od:	


mx := max(op(T));

for i from 1 to n do
	if mx = T[i] then
		return([i,T]);
	fi:
od:

ERROR(NotExpect);	
end:


# Input: positive integer n, 
# starting list of permutation S,
# number of shuffling k
# Output: the next most 
# likely elements after S
# Try: 
# FastBG(10,2,[5]); 
# SlowBG(10,2,[5]);
FastBG := proc(n,k,S) option remember;
FastBestGuess(n,2^k,S);
end:



################################
# Section 3.4: 
# Generating functions
# of n cards with k shuffles.
#
# Functions:
# CorGuess(A,k), GenK(n,k,q),
#
################################

# Input: List of permutation of 
# length n, positive integer k
# Output: number of correct guess
# under k-shuffling of deck A 
# Try: 
# CorGuess([1,2,4,3],2);
CorGuess := proc(A,k) option remember;
local i,n,count;

count := 0;
n := nops(A);

for i from 1 to nops(A) do
	if A[i] = FastBG(n,k,[op(1..i-1,A)])[1] then
		count := count+1;
	fi:	
od:
count;
end:


# Input: positive integers n,k 
# and symbolic q
# Output: Generating function 
# for k shuffles
# Try:
# GenK(3,2,q);
# GenD(11,q); GenK(11,1,q);
GenK := proc(n,k,q) option remember;
local i,A,S,D;

D := 0;
S := InSeq(n,2^k); 

for i from 1 to nops(S) do 
for A in S[i] do
	D := D+binomial(n+2^k-i,n)*q^CorGuess(A,k);
od: od:
D;
end:


################################
# Section 3.5: Eulerian numbers
#
# Functions:
# Euler1(n,m), ForEuler1(n,m), 
# IdEu(x,n),
#
#################################

# Input: Positive integer n,
# non-negative integer k
# Output: Eulerian numbers, i.e.
# number of permutation of
# length n with k descendants
# Try: 
# seq(print([seq(Euler1(n,k),k=0..n-1)]),n=1..7);
Euler1 := proc(n,k) option remember;

if k = 0 then
	return(1);
elif k > n or k < 0 then
	return(0);
fi:

(n-k)*Euler1(n-1,k-1)+(k+1)*Euler1(n-1,k);
end:
 

# Test formula for Eulerian numbers
# Input: Positive integer n,
# non-negative integer k
# Output: Eulerian numbers
# Try: 
# [seq(Euler1(7,k),k=0..6)];
# [seq(ForEuler1(7,k),k=0..6)];
ForEuler1 := proc(n,k) option remember;
local i;

add( (-1)^i*binomial(n+1,i)*(k+1-i)^n ,i=0..k);
end:


# Test Worpitzky's identity
# Input: positive integers x,n
# Output: right side and
# the difference between left and right
# side in Worpitzky's identity
# Try: IdEu(13,11);
IdEu := proc(x,n) option remember;
local k,A;

A := add( Euler1(n,k)*binomial(n+x-k-1,n)
,k=0..x-1);

[A,A-x^n];
end:




#############################
# Section 4.1: Application
# 
# Functions:
# Dora2(k), Dora3(k),
#
#############################

# Input: non-negative integer k
# Output: Expected number of # of
# correct guesses after k-shuffles
# of deck of 2 cards
# Try: [seq(Dora2(k),k=0..5)];
Dora2 := proc(k) option remember;
local m,M,P;

M := m-> binomial(2+2^k-m,2)/2^(2*k);
P := [M(1),M(2)];

2*P[1]+P[2];
end:


# Input: non-negative integer k
# Output: Expected number of # of 
# correct guesses after k-shuffles
# of deck of 3 cards
# Try: [seq(Dora3(k),k=0..5)];
Dora3 := proc(k) option remember;
local m,M,P;

M := m-> binomial(3+2^k-m,3)/2^(3*k);
P := [M(1),M(2),M(2),M(2),M(2),M(3)];

3*P[1]+2*P[2]+2*P[3]
+P[4]+2*P[5]+P[6];

end:




##################################
# Appendix A: Central binomial
#
# Functions:
# Ber(n), CenBi(k,n),
# TestCenBi(k,N),
#
##################################

# Input: non-negative integer n
# Output: Bernoulli numbers
# Try: Ber(6);
Ber := proc(n) option remember;
local t,f,B;

f := t/(exp(t)-1);
B := taylor(f,t=0,n+2);
coeff(B,t,n)*n!;
end:


# Output: Centeral binomial 
# coefficients up to the first 
# K terms computed by Bernoulli 
# numbers and Taylor series
# Input: number of terms K 
# and symbolic n
# Try: CenBi(5,n);
CenBi := proc(K,n) option remember;
local i,x,A,F;

A:= add(Ber(2*i)*x^(2*i-1)/(i*(2*i-1))
*(1-1/4^i),i=1..K/2+1);
F := taylor(exp(-A),x=0,K+2);

4^n/sqrt(Pi*n)*add( coeff(F,x,i)/n^i,i=0..K);
end:


# Test numerical values of 
# approximation of central
# binomial
# Pretty good!
# Input: number of terms K 
# and some numeric n
# Output: relative difference
# between the actual and approx.
# Try: TestCenBi(5,40.5);
TestCenBi := proc(K,n) option remember;
local N,A,B;

Digits:= 50:
A := evalf(subs(N=n,CenBi(K,N)));
B := evalf(binomial(2*n,n));
#print(A,B);
(A-B)/B;
end: