# restart; 
# read `\c:/Users/Asus/Google Drive/MyNote/ProbBridge.txt`;

with(combinat):

# Check the probability distribution
# of things I am interested in the
# game of Bridge


###################################
# Section 1: Distribution of cards
#
# Functions:
# PrCard(a,b,c,d), RunPrCard(),
#
# ParRow(n,r), Rev(P), Rep(P),
# Sort1(SS),
#
###################################

# Try: PrCard(4,3,3,3);
# PrCard(5,3,3,2);
PrCard := proc(a,b,c,d) option remember;
local A;

if a+b+c+d <> 13 then
	ERROR("BadInput");
fi:	

A := binomial(13,a)*binomial(13,b)
*binomial(13,c)*binomial(13,d)/
binomial(52,13);

evalf(A);
end:

# Try: RunPrCard();
RunPrCard := proc() option remember;
local i,P,Pro,SS;

SS := [];

for P in ParRow(13,4) do
	P := [op(P),0$(4-nops(P))];
	Pro := PrCard(op(P))*4!/Rep(P);
	SS := [op(SS),[P,Pro]];
od:

SS := Sort1(SS);
for i from 1 to nops(SS) do 
	if SS[i][2] > 0.01 then
		#print(SS[i][1],"with prob =",SS[i][2]);
	fi:	
od:	

#return( add(S[2],S in SS));
return(SS):
end:

#############################################
# Auxillary Functions

# Output: partition of n 
# with at most r parts
# Try: ParRow(6,3);
ParRow := proc(n,r) option remember;
local t,SS;

SS :=[];

for t in partition(n) do 
	if nops(t) <= r then
		SS := [op(SS),Rev(t)];
	fi:	
od:

SS;
end:

# Try: Rev([2,3,5,1]);
Rev := proc(P) option remember;
local n:

n := nops(P);
[seq(P[n-i],i=0..n-1)];
end:


# Repeat element with factorial in P
# Try: Rep([3,3,2,1,1,1]);
Rep := proc(P) option remember;
local i,m,count;

m := 1;
count := 1;

for i from 2 to nops(P) do
	if P[i] <> P[i-1] then
		count := 0;
	fi:
	count := count+1;
	m := m*count;
od:
m;
end:

# Sort by the biggest prob.
Sort1 := proc(SS) option remember;
local i,j,TT;
TT := [SS[1]];

for i from 2 to nops(SS) do
	j := 1;
	while j <= nops(TT) and SS[i][2] <= TT[j][2]  do
		j := j+1;
	od:

	TT := [op(1..j-1,TT),SS[i],op(j..nops(TT),TT)];	
od:

TT;
end:




###########################
# Section 2: Distribution 
# of my opp.'s cards
#
# Functions:
# HoldCards(N),
#
###########################
# Try: HoldCards(7);
HoldCards := proc(N) option remember;
local i,P,A,Pr;

if N<7 or N>12 then
	ERROR("NoCase");
fi:	

P := 13-N;
A := add(binomial(13,P-i)*binomial(13,i),i=0..P);

for i from 0 to P/2 do
	Pr := binomial(13,P-i)*binomial(13,i)/A;
	if i <> P-i then
		print([P-i,i], "with Prob", 2*evalf(Pr));
	else
		print([P-i,i], "with Prob", evalf(Pr));
	fi:	
od:

return();
end:

###########################
# Section 3: Simulation 
#
# Functions:
# Run1(N), Run2(n,N),
# SimSuit(),
#
###########################
# Double check the result 
# in section 1 and 2.

# Test the distribution of 
# suits in the player's hand
# Try: Run1(200000);
# RunPrCard();
Run1 := proc(N)
local i,S,P;

for P in ParRow(13,4) do
	P := [op(P),0$(4-nops(P))];
	S(P):=0;
od:

for i from 1 to N do
	P := Rev(sort(SimSuit()));
	S(P):=S(P)+1; 
od:

for P in ParRow(13,4) do
	P := [op(P),0$(4-nops(P))];
	if evalf(S(P)/N)>0.01 then
		print(P,"with probability",evalf(S(P)/N));
	fi:	
od:

return():
end:

# Test prob. dist of the opp.
# Input: number of simul.
# Output: prob. dist. when 
# team has 7 cards and up
# Try: Run2(200000);
# HoldCards(7);
Run2 := proc(N)
local i,a,P,total,S,m;

total :=[0$7];
S := [[0$4]$7];

for i from 1 to N do
	P := SimSuit():
	if P[1]+P[3] >= 7 then
		m := P[1]+P[3]-6;
		a := min(P[2],P[4]);
	else	
		m := P[2]+P[4]-6;
		a := min(P[1],P[3]);
	fi:
	total[m] := total[m]+1;
	S[m,a+1] := S[m,a+1]+1;
od:

for m from 1 to 7 do
if total[m] > 0 then
	print("For team with toal", m+6, "cards.");
	for a from 0 to (7-m)/2 do
		print([7-m-a,a], "with prob.", evalf(S[m,a+1]/total[m]));
	od:
fi:
od:

return(total,add(total[i],i=1..7));
end:


############################################
# Try: SimSuit();
SimSuit := proc() 
local roll,s,S,P;

P :=[0$4];
S := {};
roll := rand(1..52); 

while nops(S) < 13 do
	S := S union{roll()};
od:
	
for s in S do
	P[s mod 4+1] := P[s mod 4+1] +1; 
od:	

P;
end: