# restart; read `c:/Users/Asus/Google Drive
# /Aek/MomentReview/SchurT.txt`;

with(gfun, guessgf):

#print(`Welcome to SchurT, a maple package`):
#print(`written by Thotsaporn (Aek) Thanatipanonda`):
#print(`First Version: November 16, 2013 `):
#print(`Second Version: February 29, 2020 `):
#print(`Find moment of Schur triple x+y=z`):
#print(`on [1,n] interval with r-coloring.`):

Help := proc()
print(`Section1: Direct Computation`);
print(`The functions are`):
print(`Schur(n), IntSchur(S,n), LastSchur(S),`);
print(`MoSchur(n,m), MoMeanSchur(n,m).`);
print(``);
print(`Section2: Overlapping Stage Approach`);
print(`The functions are`):
print(`OSA(n,r,m), Dist(n,m), Overlap(n,m),`);
print(`Triples(n,m), Chain(A)`);
print(``);
print(`Section3: Formula`);
print(`The functions are`):
print(`ForCol1(n,m), ForMo1(n,r),`);
print(`ForMo2(n,r), ForMoMean2(n,r),`);
print(`GetForMo2(n,r),`);
print(`Expand1(SS,n), Exp(T,nn,n)`);
print(``);
print(`Section4: Guessing Formula`);
print(`The functions are`):
print(`MyGuess(N,m,n), Data1(N,m),`);
print(`GuessPeriodPol(L,x,maxDegree),`);
print(`GuessPol(L,x,maxDegree)`);
print(``);
end:



#################################
# Section 1: Direct Computation

# Input: the length of interval
# Output: distribution of number of 
# Schur triples on all of the 
# possible 2-colored interval.
# Try: Schur(5);
Schur := proc(n) option remember;
local i,t,b,T,N;

T := IntSchur([],n);
b := n*(n-1)/2;
N := array(0..b);

for i from 0 to b do
	N[i] := 0;
od:

for t in T do
	N[t] := N[t]+1;
od:

[seq(N[i],i=0..b)];
end:


# Input: the interval S and the 
# aim length n.
# Output: distribution of number of 
# Schur triples of all possible 
# 2-colored intervals of length n 
# starting with S.
# Try: IntSchur([1,1,2],5);
IntSchur := proc(S,n) option remember;
local t,s;

t := LastSchur(S);

if nops(S) = n then
	return([t]);
fi:

[seq(t+s, s in IntSchur([op(S),1],n)),
seq(t+s, s in IntSchur([op(S),2],n))];
end:


# Input: Interval S, each entry in {1,2}
# Output: number of Schur triples 
# involves the last entry.
# Try: LastSchur([1,2,1,2,2,1,1,2]);
LastSchur := proc(S) option remember;
local i,n,t,count;

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

if n <> 0 then
	t := S[n];
fi:

for i from 1 to n/2 do
	if t=S[i] and t=S[n-i]  then
		count := count+1;
	fi:
od:
count;
end:


# Average on sample space of length n
# of Schur^r by brute force.
# Try: MoSchur(5,1);
# seq(MoSchur(n,2),n=0..10);
MoSchur := proc(n,m) option remember;
local i,s,S;

S := Schur(n);
add(S[i]*(i-1)^m,i=1..nops(S))
/add(s, s in S);
end:


# Try: MoMeanSchur(5,1);
# [seq(MoMeanSchur(5,m),m=0..4)];
MoMeanSchur := proc(n,m) option remember;
local i,s,Mean,S;

S := Schur(n);
Mean := MoSchur(n,1);

add(S[i]*(i-1-Mean)^m,i=1..nops(S))
/add(s, s in S);
end:


#####################################
# Section 2: Overlapping Stage Approach

# moment of Schur Triples.
# Input: length of interval n, 
# with r coloring and m-th moment.
# Try: OSA(3,2,1);
# [seq(OSA(n,3,1),n=0..10)];
# [seq(ForMo1(n,3),n=0..10)];
# [seq(OSA(n,2,2),n=0..7)];
OSA := proc(n,r,m) option remember;
local c,s,TT;

TT := Dist(n,m);
add(r^(n-s)
*add(TT[s][c]*r^c,c=1..m)
,s=1..nops(TT))/r^n;
end:


# Distribution of Schur Triples.
# Input: length of interval n 
# and m-th moment.
# Try: Dist(3,1);
# [seq(Dist(n,2),n=0..10)];
Dist := proc(n,m) option remember;
local a,s,t,TT;

TT := [[0$m]$(3*m)];

for a in Overlap(n,m) do
	s := nops(`union`(seq(t, t in a))); 
	TT[s][Chain(a)] := TT[s][Chain(a)]+1;
od:

TT;
end:


# Try: Overlap(5,2);
Overlap := proc(n,m) option remember;
local a,s,t;

if m =1 then
	return([seq([a], a in Triples(n))]);
fi:

[seq(seq(  [op(s),op(t)]
,s in Overlap(n,1)),t in Overlap(n,m-1))];
end:


# Set of all triples on interval of length n.
# Try: Triples(3);
Triples := proc(n) option remember;
local i;

if n=0 then
	return({});
fi:

{seq({i,n-i,n},i=1..n/2), op(Triples(n-1))};
end:


# Figure out how many chian do A have?
# Try: Chain([{1,4,5},{1,4,5},{1,4,5}]);
Chain := proc(A) option remember;
local a,i,temp,T,C;

C := {};

for a in A do
	T := {};
	for i from 1 to nops(C) do 
		if a intersect C[i] <> {} then
			T := T union {i};
		fi:		
	od:
	temp := `union`(a,seq(C[i], i in T));
	C := (C minus `union`(seq({C[i]}, i in T)))
		union {temp};
od:

nops(C);
end:


#################################
# Section3: Formulas

# Formula for the m-th moment 
# of length n with 1 color.
# Try:  [seq(ForCol1(n,2),n=0..10)];
ForCol1 := proc(n,m) option remember;
	(floor(n/2)*ceil(n/2))^m;
end:


# Formula for the second moment 
# of length n with r colors.
# Try:  [seq(ForMo1(n,2),n=0..10)];
# [seq(MoSchur(n,1),n=0..10)];
ForMo1 := proc(n,r) option remember;

if n mod 2 = 1 then	
	(n-1)*(n-1+2*r)/(4*r^2);
else
	n*(n-2+2*r)/(4*r^2);
fi:

end:


# Formula for the second moment 
# of length n with r colors.
# Try: [seq(ForMo2(n,2),n=0..15)];
# [seq(MoSchur(n,2),n=0..15)];
ForMo2 := proc(n,r) option remember;
local S;

S := [ 1/48/r^4*(n-1)
	  *(3*n^3-9*n^2+12*n^2*r+24*n*r^2
	  -27*n+24*r^3-16*r^2-76*r+65),
		1/48/r^4*
	(24*n*r^3+24*n^2*r^2-16*n*r^2-64*r^2
	-112*n*r+128*r-12*n^3-24*n^2+104*n
	-64+12*n^3*r+3*n^4),
		1/48*(24*n*r^3-24*r^3+48*r^2-40*r^2*n
	+24*r^2*n^2-12*r*n^2-76*r*n+12*r
	-12*n^3-18*n^2+92*n-33+12*r*n^3+3*n^4)/r^4,
		1/48*(24*n*r^3-32*r^2-16*r^2*n
	+24*r^2*n^2+64*r-112*r*n-12*n^3
	-24*n^2+104*n-32+12*r*n^3+3*n^4)/r^4,
		1/48*(24*n*r^3-24*r^3-16*r^2
	-40*r^2*n+24*r^2*n^2+140*r-76*r*n
	-12*r*n^2-97+92*n-18*n^2-12*n^3
	+12*r*n^3+3*n^4)/r^4,
		1/48*n*(24*r^3-16*r^2+24*r^2*n
	-112*r+104-24*n-12*n^2+12*r*n^2
	+3*n^3)/r^4,
		1/48*(n-1)*(24*r^3+24*r^2*n
		-16*r^2-76*r-9*n^2-27*n+65
		+12*r*n^2+3*n^3)/r^4,
		1/48*(24*n*r^3+24*r^2*n^2-16*r^2*n
	-64*r^2-112*r*n+128*r-12*n^3
	-24*n^2+104*n-64+12*r*n^3+3*n^4)/r^4,
		1/48*(24*n*r^3-24*r^3+48*r^2
	-40*r^2*n+24*r^2*n^2-12*r*n^2
	-76*r*n+12*r-12*n^3-18*n^2+92*n
	-33+12*r*n^3+3*n^4)/r^4,
		1/48*(24*n*r^3-32*r^2-16*r^2*n
	+24*r^2*n^2+64*r-112*r*n-12*n^3
	-24*n^2+104*n-32+12*r*n^3+3*n^4)
	/r^4,
		1/48*(24*n*r^3-24*r^3-16*r^2-40*r^2*n
	+24*r^2*n^2+140*r-76*r*n-12*r*n^2-97
	+92*n-18*n^2-12*n^3+12*r*n^3+3*n^4)
	/r^4,
		1/48*n*(24*r^3-16*r^2+24*r^2*n-112*r
	+104-24*n-12*n^2+12*r*n^2+3*n^3)/r^4];

if n mod 12 = 0 then	
	S[12];
else
	S[n mod 12];
fi:

end:


# Formula for the second moment 
# about mean of length n with r colors.
# Try: [seq(ForMoMean2(n,2),n=0..10)];
ForMoMean2 := proc(n,r) option remember;
local S;

S := [	1/12*(n-1)*
(-9*n+6*n*r+3*n*r^2+17-22*r-r^2+6*r^3)/r^4, 
		1/12*(-9*n^2+6*n^2*r+3*n^2*r^2+6*n*r^3
-4*n*r^2-16*r^2-28*n*r+32*r+26*n-16)/r^4, 
		1/12*(6*r+26*n+9*r^2-4*n*r^2+6*n*r^3
	+3*n^2*r^2-28*n*r-6*r^3+6*n^2*r-9*n^2-9)/r^4, 
		1/12*(-9*n^2+6*n^2*r+3*n^2*r^2+6*n*r^3
	-8*r^2-4*n*r^2+16*r-28*n*r+26*n-8)/r^4, 
		1/12*(38*r+26*n-7*r^2-4*n*r^2+6*n*r^3
	+3*n^2*r^2-28*n*r-6*r^3+6*n^2*r-9*n^2-25)/r^4, 
		1/12*n*(-9*n+6*n*r+3*n*r^2+6*r^3-4*r^2
	-28*r+26)/r^4, 
		1/12*(n-1)*(-9*n+6*n*r+3*n*r^2+17
	-22*r-r^2+6*r^3)/r^4, 
		1/12*(-9*n^2+6*n^2*r+3*n^2*r^2
	+6*n*r^3-4*n*r^2-16*r^2-28*n*r+32*r
	+26*n-16)/r^4, 
		1/12*(6*r+26*n+9*r^2-4*n*r^2+6*n*r^3
	+3*n^2*r^2-28*n*r-6*r^3+6*n^2*r-9*n^2-9)/r^4, 
		1/12*(-9*n^2+6*n^2*r+3*n^2*r^2+6*n*r^3
	-8*r^2-4*n*r^2+16*r-28*n*r+26*n-8)/r^4, 
		1/12*(38*r+26*n-7*r^2-4*n*r^2+6*n*r^3
	+3*n^2*r^2-28*n*r-6*r^3+6*n^2*r-9*n^2-25)/r^4, 
		1/12*n*(-9*n+6*n*r+3*n*r^2+6*r^3
	-4*r^2-28*r+26)/r^4];

if n mod 12 = 0 then	
	S[12];
else
	S[n mod 12];
fi:

end:



###############
# Sub-function
###############
# Formulas for the second moment 
# of length n with r colors.
# Try: GetForMo2(n,r);
GetForMo2 := proc(n,r) option remember;
local s,i,j,k,S;

S := Expand1(MyGuess(72,2,n),n);
k := nops(S[1][1]);

[seq(factor(simplify(
add(add( S[s][i][j]*r^(i-s)
,i=1..2),s=1..6))),j=1..k)];
end:


# Try: Expand1([[[n-1/2,n/2],[n]],[[0$3],[0$2]]],n)
Expand1 := proc(SS,n)
local i,j,s1,s2,nn;

s1 := nops(SS);
s2 := nops(SS[1]);
nn := lcm(seq(seq(nops(SS[i][j]),j=1..s2),i=1..s1));

[seq([seq( Exp1(SS[i][j],nn,n)
,j=1..s2)],i=1..s1)];
end:


# Try: Exp([(n-1)/2,n/2],6,n);
Exp1 := proc(T,nn,n) option remember;
local i;

[seq(op(T),i=1..nn/nops(T))];
end:


##############################
# Section 4: Guessing formula

# Try: MyGuess(72,2,n);
MyGuess := proc(N,m,n) option remember;
local DD;

DD := Data1(N,m);

[seq(
[GuessPeriodPol(DD[i][1],n,m*3),
GuessPeriodPol(DD[i][2],n,m*3)]
,i=1..nops(DD))];

end:


# I can do guessgf to get gf from these datas.
# Try: Data1(30,2);
Data1 := proc(N,m) option remember;
local n,SS,T;

SS := [[[]$m]$(3*m)];

for n from 1 to N do 
	T := Dist(n,m);
	SS := [seq( 
		  [[op(SS[i][1]),T[i][1]],
		  [op(SS[i][2]),T[i][2]]]
		,i=1..3*m)];
od:

SS;
end:


###########
# Guessing
###########
# GuessPeriodPol(L,x, maxDegree) 
# - Find Polynomial of Period of list 
# Input: a list of numbers
# Output: a list polynomial P of x, 
# of degree min(nops(L)-2,max)
# such that L[i]=P(i), for i=1..nops(L) 

# Try: LL:=[0,0,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6]; 
# GuessPeriodPol(LL,n,3);
GuessPeriodPol := proc(list, x, maxDegree)
local i, j, found, period, interval, result, ret;

found := false;

for period from 1 to nops(list)/4 while not(found) do
	ret   := [];
	found := true;

	for i from 1 to period while found do
		interval := [seq(list[j*period+i]
		, j=0..(nops(list)-i)/period)];
		result   := factor(subs(x=(x-i)/period+1
		, GuessPol(interval, x, maxDegree)));
		ret      := [op(ret), result];
		found    := (result <> FAIL);
	od;
od;

if found then
	return ret;
else 
	return FAIL;
fi;
end:


# GuessPol(L,x, maxDegree) - Guess Polynomial 
# Input: a list of numbers, vairable x and maxDegree.
# Output: a polynomial P of x, 
# of degree min(nops(L)-2,max)
# such that L[i]=P(i), for i=1..nops(L)

# Try: GuessPol([1,3,6,10,15,21,28,36],n,2);
GuessPol := proc(L, x, maxDegree) 
local d, P, a, var, eq;
	
d := min(nops(L)-2, maxDegree);
P := add(a[i]*x^i, i=0..d);
var := {seq(a[i], i=0..d)};
eq := {seq(subs(x=i, P)=L[i], i=1..nops(L))};
var := solve(eq, var);

if var = NULL then
	return FAIL;
else
	return factor(subs(var, P));
fi;
end: