#its ok to post
#TaerimKim,10/23/2020,Assignment 12

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

#1. this seems relatively easy to construct 
#its just sequencing all L[i] elements in form of CP(L[i],L[j}) using the same logic as CP

CPg:=proc(L) local i, j, l, Lf;
l:=nops(L);
Lf:=[];

for i from 1 to l-1 do
for j from i+1 to l do
Lf:=[op(Lf),CP(L[i],L[j])];
od;
od;

Lf;
end;

#lets check

CP([45, 65, 47], [43, 43, 65]);
           [88, 88, 110, 108, 108, 130, 90, 90, 112]
#this have to be same as CPg([[45, 65, 47], [43, 43, 65]])

CPg([[45, 65, 47], [43, 43, 65]])
[[88, 88, 110, 108, 108, 130, 90, 90, 112]]

#true

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

#2. the difference between AveClever/kthMomentClever and new formula is that we don't need to incorportate WtEn in the logic

#so, by simplifying the formula by excluding WtEn gives us

AveGF:=proc(f,x); 
subs(x=1,diff(f,x))/subs(x=1,f):
end:

kthMomentGF:=proc(f,x,k) local i,mu,f1; 
mu:=AveGF(f,x):
f1:=f/x^mu/subs(x=1,f):
for i from 1 to k do
f1:=expand(x*diff(f1,x)):
od:
subs(x=1,f1):
end:

#Let's see if that is true

evalb(AveGF(WtEn(L, x), x) = AveClever(L));
                              true

evalb(kthMomentGF(WtEn(L,x),x,k)= kthMomentClever(L,k)) #here k needs to be numeric so we are inputing k=2
evalb(kthMomentGF(WtEn(L, x), x, 2) = kthMomentClever(L, 2));
                              true


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

#3. 
#we know that kthmoment of the mean can be acquired by our new formula kthMomentGf with f,x as input

ScaledMomentGF:=proc(f,x,k);
kthMomentGF(f,x,k)/(kthMomentGF(f,x,2)^(k/2));
end;

#to check, Lets use the L as [45, 65, 47] and Weight-enumerator is 
WtEn([45, 65, 47], x);
                        x^65 + x^47 + x^45
f:=%

ScaledMomentGF(f,x,6)
8179931/3014284

#compare that with original formula

kthMomentGF(f,x,6)/(kthMomentGF(f,x,2)^(6/2))
8179931/3014284

#same, so its good

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

#4. lets formulate the given 

fu:=proc(n) ; 
seq(ScaledMomentGF(((1+x)/2)^n,x,i),i=2..10);
end;

#the explicit form of this sequence should be fu with n as an input
fu(n)
1,0,(16 (-1/8 n+3/16 n^2))/(n^2),0,(64 (1/4 n+15/64 n^3-15/32 n^2))/(n^3),0,(256 (-17/16 n+105/256 n^4-105/64 n^3+147/64 n^2))/(n^4),0,1024*(31/4*n + 945/1024*n^5 - 1575/256*n^4 + 4095/256*n^3 - 1185/64*n^2)/n^5
u:=[%] #going to put the sequence into a set

#then, going to sequence this by testing limit to infinity to every element in the set
seq(limit(u[i], n = infinity), i = 1 .. 9)
1, 0, 3, 0, 15, 0, 105, 0, 945

#lets look this up in OEIS
#the set is called A123023 and is the sequence of a(n) = (n-1)*a(n-2), a(0)=1, a(1)=0, which is number of ways of separating n terms into pairs.