# restart; read `c:/Users/Asus/Google Drive/Aek/
# AverageGame/SimpDre.txt`;

with(SolveTools):


##################################
# Section 1: Simplify Dreidel:
# one to the pot or Take the whole pot
#
##################################
# Section 1.1: Solve values directly
# 
# Functions:
# Sim2(), Sim3(),
# SDSolve(p1,pot,p2), MakeVar(A,T),
# SD(p1,pot,p2), 
#
##################################

# with N=2
# Try: Sim2();
Sim2 := proc() option remember;
local T,eq,var;

eq := { 
T[1,2,1] = 1+1/2*(T[1,3,0]+T[0,2,2]),
T[1,3,0] = 0,
T[0,2,2] = 0
};
var := {T[1,2,1],T[1,3,0],T[0,2,2]};
Linear(eq,var);
end:

# Variant: one to the pot 
# or Take the whole pot
# with N=3 i.e. 3 tokens each 
# at the start
# Try: Sim3();
Sim3 := proc() option remember;
local T,eq,var;

eq := { 
T[2,2,2] = 1+1/2*(T[2,3,1]+T[1,2,3]),
T[2,3,1] = 1+1/2*(T[1,4,1]+T[0,2,4]),
T[1,2,3] = 1+1/2*(T[3,3,0]+T[2,2,2]),
T[1,4,1] = 1+1/2*(T[1,5,0]+T[0,2,4]),
T[0,2,4] = 0,
T[3,3,0] = 0,
T[1,5,0] = 0
};
var := {T[2,2,2],T[2,3,1],T[1,2,3],T[1,4,1],
T[0,2,4],T[3,3,0],T[1,5,0]};

Linear(eq,var);
end:



# Variant with general starting tokens and pot
# Try: SDSolve(2,2,2);
# [seq( subs(SDSolve(n-1,2,n-1),T[n-1,2,n-1]),n=2..16)];
SDSolve := proc(p1,pot,p2) option remember;
local i,s,A,SS,NS,eq,var;

eq := { };
var := {};
SS := {T[p1,pot,p2]};

while SS <> {} do 
	NS := {};
	for s in SS do
		if [op(s)][1]=0 or [op(s)][3]=0 then
			eq := eq union {s=0};
		else	
			A := MakeVar(s,T); 
			eq := eq union {s=1+1/2*add(A[i],i=1..2)};
			NS := NS union {seq(A[i],i=1..2)};
		fi:
		
		var := var union {s};
	od:	
	SS := NS minus var; 
od:

Linear(eq,var);
end:


# Make choice for "SDSolve" problem
# Input: Current state of SDSolve problem
# Output: the possible state: 
# pay 1 or take all

# Try: MakeVar(T[2,1,3],T);
MakeVar := proc(A,T) option remember;
local p1,pot,p2;

p1  := [op(A)][1];
pot := [op(A)][2];
p2  := [op(A)][3];

[T[p2,pot+1,p1-1],T[p2-1,2,p1+pot-1]];
end:

# Value of the actual points
# Try: 
# N:=10: [seq(SD(i,2,N-i),i=0..N)];
# A:=[seq(SD(1,p,3),p=2..10)];
SD := proc(p1,pot,p2) option remember;
subs(SDSolve(p1,pot,p2),T[p1,pot,p2]);
end:


####################################
# Section 1.2: 
# Another idea on C and CP
# and solve for T(a,b) directly
#
# Functions:
# GDela2(a,p,b), GDelb2(a,p,b),
# GenT(a,b,T), AllT(n), 
# MyT(a,b),
#
# MyCP(a,p,b),
#
####################################
# Toward solution

# T(a+2,b)-2*T(a+1,b)+T(a,b)
# Go to 0 as a,b big
# Try: evalf([seq(GDela2(a,2,13),a=1..20)]);
GDela2 := proc(a,p,b) option remember;
SD(a+2,p,b)-2*SD(a+1,p,b)+SD(a,p,b);
end:


# T(a,b+2)-2*T(a,b+1)+T(a,b)
# Go to 0 as a,b big
# Try: evalf([seq(GDelb2(a,2,3),a=1..20)]);
GDelb2 := proc(a,p,b) option remember;
SD(a,p,b+2)-2*SD(a,p,b+1)+SD(a,p,b);
end:


# Solving Del directly using SD
# Generate T(a,b) with pot =2
# Try:  GenT(1,4,T);
# [seq(GenT(a,4,T)-SD(a,2,4),a=1..10)];
GenT := proc(a,b,T) option remember;
local i;

[add(1/2^(2*i-1)*T(b-i,a+i),i=1..min(a,b))
,add(1/2^(2*i-2)*T(a-i,b+i),i=2..min(a,b+1))
,add(1/2^i,i=0..min(2*(a-1),2*b-1))];

end:


# Faster version to solve the system
# Try: AllT(6);
AllT := proc(n) option remember;
local a,T,S,eq,var;

S := {seq(T(a,n-a)=T[a],a=0..n)};
eq := {seq(T[a]=add(GenT(a,n-a,T)[i],i=1..3),a=0..n)};
var := [seq(T[a],a=0..n)];

subs(Linear(subs(S,eq),var),var);
end:


# Faster way to find SD with pot=2
# Try: 
# SD(13,2,14); MyT(13,14);
MyT := proc(a,b) option remember;
AllT(a+b)[a+1];
end:

###################################################
# Another idea!!

# Shift in p with any pot
# Try: 
# evalf([seq(MyCP(a,5,45)
# -MyCP(a-1,5,45)-4/19,a=30..33)]);
# evalf([seq(MyCP(50,2,b)-4/19*50-8/19*b,b=40..43)]);
# evalf([seq(MyCP(17,11,b)-MyCP(17,2,b),b=20..27)]);
MyCP := proc(a,p,b) option remember;
SD(a,p+1,b)-SD(a,p,b);
end:



####################################
# Section 1.3: 
# Plugin equations to solve for 
# coeffs of T(a,b) 
#
# Functions:
# Case1(a,b), Case2(a,b),
#
####################################

# Case a >= b+1:
# Try: Case1(a,b);
Case1 := proc(a,b) option remember;
local T,c3,c2,c1,c0,L1,L2,L3,M,eq;

T := (a,b)->c9*a^2*b+c8*a*b^2+c7*a^3+c6*b^3
+c5*a^2+c4*b^2+c3*a*b+c2*a+c1*b+c0:
L1 := sum(1/2^i,i=0..2*b-1);
L2 := sum(T(b-i,a+i)/2^(2*i-1),i=1..b);
L3 := sum(T(a-i,b+i)/2^(2*i-2),i=2..b+1);
M:=subs({4^(-b)=0,2^(-2*b+2)=0},simplify(expand(L1+L2+L3)));

eq:=[
c9=coeff(coeff(M,b,1),a,2),
c8=coeff(coeff(M,b,2),a,1),
c7=coeff(coeff(M,b,0),a,3),
c6=coeff(coeff(M,b,3),a,0),
c5=coeff(coeff(M,b,0),a,2),
c4=coeff(coeff(M,b,2),a,0),
c3=coeff(coeff(M,b,1),a,1),
c2=coeff(coeff(M,a,1),b,0),
c1=coeff(coeff(M,b,1),a,0),
c0=coeff(coeff(M,b,0),a,0)];
solve(eq,{c9,c8,c7,c6,c5,c4,c3,c2,c1,c0});
end:


# Case a <= b:
# Try: Case2(a,b);
Case2 := proc(a,b) option remember;
local T,c3,c2,c1,c0,L1,L2,L3,M;

T := (a,b)->c3*a*b+c2*a+c1*b+c0:
L1 := sum(1/2^i,i=0..2*a-2);
L2 := sum(T(b-i,a+i)/2^(2*i-1),i=1..a);
L3 := sum(T(a-i,b+i)/2^(2*i-2),i=2..a);
M:=subs({4^(-a)=0,2^(-2*a+2)=0},simplify(expand(L1+L2+L3)));
[coeff(coeff(M,b,1),a,1),coeff(coeff(M,a,1),b,0),
coeff(coeff(M,b,1),a,0),coeff(coeff(M,b,0),a,0)];
end:



#################################
# Section 1.4:
# Least Square method for
# guessing solution.
#
# Functions:
# LSq(X,Y,x), LSSimDr(x,y), 
# SolT(a,b), SolSimDr(a,p,b),
# GuessFrac(c,N),
#
#################################

# Least square method for linear funciton
# Try: problem from chapter 8, 
# numerical analysis, Burden & Faires
# X:=[seq(i,i=1..10)]:
# Y:=[1.3,3.5,4.2,5,7,8.8,10.1,12.5,13,15.6]:
# LSq(X,Y,x);
LSq := proc(X,Y,x) option remember;
local i,m,a0,a1;

m := nops(X);
a0 := (add(X[i]^2,i=1..m)*add(Y[i],i=1..m) 
-add(X[i]*Y[i],i=1..m)*add(X[i],i=1..m))/
(m*add(X[i]^2,i=1..m)-add(X[i],i=1..m)^2);

a1 := (m*add(X[i]*Y[i],i=1..m) 
-add(X[i],i=1..m)*add(Y[i],i=1..m))/
(m*add(X[i]^2,i=1..m)-add(X[i],i=1..m)^2);

a1*x+a0;
end: 

# Try: S:=LSSimDr(a,b);
LSSimDr := proc(x,y) option remember;
local i,j,low,top,T,S,A,B1,B2;

B1 := [];
B2 := [];
low := 30:
top := 60:
T := [seq(i,i=low..top)];

for j from low to top do 	
	A := [seq(MyT(i,j),i=low..top)];
	S := evalf(LSq(T,A,x));
	if j mod 10 =0 then
		print(j,S);
	fi:	
	B1 := [op(B1),coeff(S,x,1)];
	B2 := [op(B2),coeff(S,x,0)];
od:	

expand(evalf(LSq(T,B1,y)*x+LSq(T,B2,y)));	
end:

# Solution to starting positions
# of simplied dreidel
# Try: 
# SolT(100,100); evalf(MyT(100,100));
# SolT(100,30); evalf(MyT(100,30));
SolT := proc(a,b) option remember;
local c2,c0;

c2 := -0.304636563075539048549763880073;
c0 := 2.13102617414184839703730983871;

12/19*a*b+c2*a+(c2+2/19)*b+c0;
end:


# Solution to simplied dreidel
# Try: 
# SolSimDr(15,5,30); evalf(SD(15,5,30));
SolSimDr := proc(a,p,b) option remember;
local c2,c0;

c2 := -0.304636563075539048549763880073;
c0 := 2.13102617414184839703730983871;

12/19*a*b+(c2+4*p/19-8/19)*a+(c2+8*p/19-14/19)*b
+c0+(p-2)*(c2-18/19);
end:


# Input: number c and the bound N
# Output: Guess whether the c close
# to x/n for some n<= N
# Try:
GuessFrac := proc(c,N) option remember;
local i;

for i from 1 to N do
	if c*i-floor(c*i) < 10^(-5) 
	or ceil(c*i)-c*i < 10^(-5) then
		print(i,c*i);
	fi:
od:	

return(Done);
end:



######################################
# Section 2:  Dreidel 
#
# Functions:
# Dreidel2(), Dreidel3(),
# DrSolve(p1,pot,p2), MakeCh(A,T),
# Dr(p1,pot,p2), 
# LSFullDr(x,y), SolDr(a,b);
#
#####################################


# Dreidel with N=2
# Try: Dreidel2();
Dreidel2 := proc() option remember;
local T,eq,var;

eq := { 
T[1,2,1]=1+1/4*(T[0,2,2]+T[1,1,2]+T[1,2,1]+T[1,3,0]),
T[0,2,2]=0,
T[1,3,0]=0,
T[1,1,2]=1+1/4*(T[1,2,1]+T[2,1,1]+T[2,1,1]+T[2,2,0]),
T[2,1,1]=1+1/4*(T[0,2,2]+T[1,1,2]+T[1,1,2]+T[1,2,1]),
T[2,2,0]=0
};
var := {T[1,2,1],T[0,2,2],T[1,3,0],T[1,1,2],T[2,1,1],T[2,2,0]};

Linear(eq,var);
end:

# Dreidel with N=3
# Try: Dreidel3();
Dreidel3 := proc() option remember;
local T,eq,var;

eq := { 
T[2,2,2]=1+1/4*(T[1,2,3]+T[2,1,3]+T[2,2,2]+T[2,3,1]),
T[1,2,3]=1+1/4*(T[2,2,2]+T[3,1,2]+T[3,2,1]+T[3,3,0]),
T[2,1,3]=1+1/4*(T[2,2,2]+T[3,1,2]+T[3,1,2]+T[3,2,1]),
T[2,3,1]=1+1/4*(T[0,2,4]+T[1,2,3]+T[1,3,2]+T[1,4,1]),
T[3,1,2]=1+1/4*(T[1,2,3]+T[2,1,3]+T[2,1,3]+T[2,2,2]),
T[3,2,1]=1+1/4*(T[0,2,4]+T[1,1,4]+T[1,2,3]+T[1,3,2]),
T[3,3,0]=0,
T[0,2,4]=0,
T[1,3,2]=1+1/4*(T[1,2,3]+T[2,2,2]+T[2,3,1]+T[2,4,0]),
T[1,4,1]=1+1/4*(T[0,2,4]+T[1,2,3]+T[1,4,1]+T[1,5,0]),
T[1,1,4]=1+1/4*(T[3,2,1]+T[4,1,1]+T[4,1,1]+T[4,2,0]),
T[2,4,0]=0,
T[1,5,0]=0,
T[4,1,1]=1+1/4*(T[0,2,4]+T[1,1,4]+T[1,1,4]+T[1,2,3]),
T[4,2,0]=0
};

var := {T[2,2,2],T[1,2,3],T[2,1,3],T[2,3,1],T[3,1,2],T[3,2,1],
T[3,3,0],T[0,2,4],T[1,3,2],T[1,4,1],T[1,1,4],T[2,4,0],
T[1,5,0],T[4,1,1],T[4,2,0]};

Linear(eq,var);
end:


# Solve Dreidel with a,b nuts 
# at p nuts in pots 
# Try: DrSolve(3,2,3);
DrSolve := proc(a,p,b) option remember;
local i,s,A,SS,NS,eq,var;

eq := {};
var := {};
SS := {T[a,p,b]};

while SS <> {} do 
	NS := {};
	for s in SS do
		if [op(s)][1]=0 or [op(s)][3]=0 then
			eq := eq union {s=0};
		else	
			A := MakeCh(s,T); 
			eq := eq union {s=1+1/4*add(A[i],i=1..4)};
			NS := NS union {seq(A[i],i=1..4)};
		fi:
		
		var := var union {s};
	od:	
	SS := NS minus var; 
od:

Linear(eq,var);
end:



# Make choice
# Input: Current state of dreidel
# Output: the possible 4 states
# Try: MakeCh(T[2,1,3],T);
MakeCh := proc(A,T) option remember;
local p1,pot,p2;

p1  := [op(A)][1];
pot := [op(A)][2];
p2  := [op(A)][3];

[T[p2-1,2,p1+pot-1],T[p2,ceil(pot/2),p1+floor(pot/2)]
,T[p2,pot,p1],T[p2,pot+1,p1-1]];
end:


# Give value of Dreidel at exact position
# Try: Dr(10,2,10);
# evalf([seq(Dr(a+1,2,13)-Dr(a,2,13),a=12..16)]);
# evalf([seq(Dr(13,2,b+1)-Dr(13,2,b),b=12..16)]);
Dr := proc(p1,pot,p2) option remember;
if p1 <0 or p2 < 0 or pot <0 then
	return(Error);
elif pot = 0 then
	subs(DrSolve(p1-1,2,p2-1),T[p1-1,2,p2-1]);
else
	subs(DrSolve(p1,pot,p2),T[p1,pot,p2]);	
fi:

end:


# Least squared for the full Dreidel
# Try: S:=LSFullDr(a,b);
# subs({a=35,b=22},S); evalf(Dr(35,2,22));
LSFullDr := proc(x,y) option remember;
local i,j,low,top,T,S,A,B1,B2;

B1 := [];
B2 := [];
low := 15:
top := 25:
T := [seq(i,i=low..top)];

for j from low to top do 	
	A := [seq(Dr(i,2,j),i=low..top)];
	S := evalf(LSq(T,A,x));
	if j mod 5 =0 then
		print(j,S);
	fi:	
	B1 := [op(B1),coeff(S,x,1)];
	B2 := [op(B2),coeff(S,x,0)];
od:	

expand(evalf(LSq(T,B1,y)*x+LSq(T,B2,y)));	
end:


# Approximation to dreidel
# Try: 
# SolDr(20,30); evalf(Dr(20,2,30));
# SolDr(10,10)*10/60;
# SolDr(20,20)*10/60;
SolDr := proc(a,b) option remember;
local c3,c2,c1,c0;

c3:= 2.21814151862618181904832628843;
c2 := -1.09709667033405910669478639669;
c1 := -0.447079544643588135688652268182;
c0 := 2.83880783734231869675987135868;

c3*a*b+c2*a+c1*b+c0;
end:




##########################
# Section 2.1: Testing!!!
#
# Functions: 
# ForDr1(a,b), ForPot(a,p,b), 
# Help1(a,p,b), 
#
# CoEq(p,c1,c2,c3,c4), 
# MyEq(a,p,b), Solve2(),
#
##########################

# Formula for dreidel when pot =2
# Try:
# [seq(ForDr1(15,b)-Dr(15,2,b),b=13..16)];
ForDr1 := proc(a,b) option remember;
local i;

+add(1/4^i, i=0..min(2*a-2,2*b-1))
+add(Dr(b-i,2,a+i)/4^(2*i-1),i=1..min(a,b))
+add(Dr(a-i,2,b+i)/4^(2*i-2),i=2..min(a,b+1))
+add(Dr(b-i,i+1,a+1)/4^(2*i+1),i=0..min(a,b)-1)
+add(Dr(a-i,i+1,b+1)/4^(2*i),i=1..min(a,b+1)-1)
+add(Dr(b-i,2+2*i,a-i)/4^(2*i+1),i=0..min(a,b)-1)
+add(Dr(a-i,1+2*i,b+1-i)/4^(2*i),i=1..min(a,b+1)-1)
;
end:


# Formula for dreidel for general pot
# Try:
# [seq(ForPot(15,4,b)-Dr(15,4,b),b=13..16)];
ForPot := proc(a,p,b) option remember;
local i;

+add(1/4^i, i=0..min(2*a-2,2*b-1))
+add(Dr(b-i,2,a+i+p-2)/4^(2*i-1),i=1..min(a,b))
+add(Dr(a-i,2,b+i+p-2)/4^(2*i-2),i=2..min(a,b+1))
+add(Dr(b-i,i+ceil(p/2),a+floor(p/2))/4^(2*i+1)
,i=0..min(a,b)-1)
+add(Dr(a-i,i+floor(p/2),b+ceil(p/2))/4^(2*i)
,i=1..min(a,b+1)-1)
+add(Dr(b-i,p+2*i,a-i)/4^(2*i+1),i=0..min(a,b)-1)
+add(Dr(a-i,p-1+2*i,b+1-i)/4^(2*i),i=1..min(a,b+1)-1)
;

end:


# Generate the recurrence
# Take whole pot, take smaller half,
# stay and pay 1 to the pot
# Try: Help1(11,2,3);
Help1 := proc(a,p,b) option remember;

if a = 0 or b =0 then
	return();
fi:	

print(  [b-1,2,a+p-1],[b,ceil(p/2),a+floor(p/2)]
,[b,p,a], [b,p+1,a-1] );
Help1(b,p+1,a-1);
end:





################################################

# Output: List of coeffs of eq. 
# when pot=p
# Try: CoEq(3,c1,c2,c3,c4); 
CoEq := proc(p,c1,c2,c3,c4) option remember;
local Tp,Th;

if p=1 then
	[c1,-c1/2+c2/3+2*c3/3,-c1/2+2*c2/3+c3/3,2-c2/2-c3/2+c4];
elif p=2 then
	[c1,c2,c3,c4];
else
	Tp := CoEq(p-1,c1,c2,c3,c4);
	Th := CoEq(ceil((p-1)/2),c1,c2,c3,c4);
	[c1,
	3*c1-c1*p-c1*floor((p-1)/2)+4*Tp[3]-Th[2]-Tp[2]-c2, 
	4*Tp[2]+c1-c3-Th[3]-Tp[3], 
	-4+4*Tp[2]+c3+c2-c1-c4-Tp[3]+3*Tp[4]
	-Th[3]-Th[4]+c1*p-c3*p-Th[3]*floor((p-1)/2)];
fi:

end:


# Settle eq for each p
# Try:
# [seq(MyEq(a,1,25)-Dr(a,1,25),a=21..23)];
# [seq(MyEq(a,3,25)-Dr(a,3,25),a=21..23)];
# [seq(MyEq(a,15,25)-Dr(a,15,25),a=11..13)];
MyEq := proc(a,p,b) option remember;
local s1,s2,s3,s4,c1,c2,c3,c4,A;

Digits:=30:
c1 := 2.21814151862618181904832628843;
c2 := -1.09709667033405910669478639669;
c3 := -0.447079544643588135688652268182;
c4 := 2.83880783734231869675987135868;

A := CoEq(p,s1,s2,s3,s4);
subs({s1=c1,s2=c2,s3=c3,s4=c4},
A[1]*a*b+A[2]*a+A[3]*b+A[4]);
	
end:



# Find constants from relation
# between 1 and 2 in the pot 
# Try: Solve2();
Solve2 := proc() option remember;
local a,b,T2,T1,M;

T2 := (a,b)->c1*a*b+c2*a+c3*b+c4:
T1 := (a,b)->d1*a*b+d2*a+d3*b+d4:

M := expand(1+1/4*(T2(b-1,a)
+T1(b,a)+T1(b,a)+T2(b,a-1)));

[d1=coeff(coeff(M,b,1),a,1),
d2=coeff(coeff(M,a,1),b,0),
d3=coeff(coeff(M,b,1),a,0),
d4=coeff(coeff(M,b,0),a,0)];
end: