#Nathan Fox
#Homework 4
#I give permission for this work to be posted online

#Read packages
with(`combinat`):

#Read procedures from class
read(`C4.txt`):

#Help procedure
Help:=proc():
    print(` ExtractData(A, i, j, a) , sdS(L) , fp(pi), MomFP(n, k) `):
    print(` HowManyGames(A, m, n) `):
end:

##PROBLEM 2##

#ExtractData(A, i, j, a): inputs a list of lists of statistical
#data, like in ClassStat() of C4.txt and outputs the sublist
#consisting only of those items for which the i-th component
#equals a, and lists only the j-th component.
ExtractData:=proc(A, i, j, a) local k, ret:
    ret:=[]:
    for k from 1 to nops(A) do
        if A[k][i] = a then
            ret:=[op(ret), A[k][j]]:
        fi:
    od:
    return ret:
end:

#Simple standard deviation
sdS:=proc(L)
    return sd([1/nops(L)$nops(L)], L):
end:

##The (simple) average height of females in our class
#AveS(ExtractData(ClassStat(), 2, 4, F)); returned 165.1666667
##The (simple) standard deviation of height of females in our class
#sdS(ExtractData(ClassStat(), 2, 4, F)); returned 6.094168435
##The (simple) average height of males in our class
#AveS(ExtractData(ClassStat(), 2, 4, M)); returned 176.3636364
##The (simple) standard deviation of height of males in our class
#sdS(ExtractData(ClassStat(), 2, 4, M)); returned 5.261744045
##The (simple) average age of females in our class
#AveS(ExtractData(ClassStat(), 2, 3, F)); returned 25.
##The (simple) standard deviation of age of females in our class
#sdS(ExtractData(ClassStat(), 2, 3, F)); returned 1.914854215
##The (simple) average age of males in our class
#AveS(ExtractData(ClassStat(), 2, 3, M)); returned 29.63636364
##The (simple) standard deviation of age of males in our class
#sdS(ExtractData(ClassStat(), 2, 3, M)); returned 12.85777251
##The (simple) average age of brown-eyed people in our class
#AveS(ExtractData(ClassStat(), 5, 3, Br)); returned 29.45454545
##The (simple) standard deviation of age of brown-eyed people in our class
#sdS(ExtractData(ClassStat(), 5, 3, Br)); returned 12.49991735
##The (simple) average age of blue-eyed people in our class
#AveS(ExtractData(ClassStat(), 5, 3, Bl)); returned 25.40000000
##The (simple) standard deviation of age of blue-eyed people in our class
#sdS(ExtractData(ClassStat(), 5, 3, Bl)); returned 5.276362384

##PROBLEM 3##

#fp(pi): number of fixed points of permutation pi
fp:=proc(pi) local i, count:
    count:=0:
    for i from 1 to nops(pi) do
        if pi[i] = i then
            count:=count+1:
        fi:
    od:
    return count:
end:

#MomFP(n, k): kth moment of the random variable, "number of fixed
#points" defined over the set of permutations on n elements
MomFP:=proc(n, k) local L, i:
    L:=permute(n):
    L:=[seq(fp(L[i]), i=1..nops(L))]:
    return MomS([1/nops(L)$nops(L)], L, k):
end:

#The first moments are all 1
#MomFP(2, 1) is 1; the other second moments are 2
#Conjecture: The expected number of fixed points is 1
#Conjecture: The variance is 2 (unless n=1)
#Conjecture (based on further experimentation): MomFP(n, k)=B(n)
#if n>=k, where B(n) is the nth Bell number

##PROBLEM 4##

#HowManyGames(A, m, n): inputs a set of positive integers, A, and
#two positive integers, m and n, and outputs the number of possible
#scoring histories that could have lead to the score m:n in a game 
#where the `atomic' scores are taken from the set A
HowManyGames:=proc(A, m, n) local f, a:
    f:=1/(1-add(x^a, a in A)-add(y^a, a in A)):
    return coeff(taylor(normal(coeff(taylor(f,x=0,m+1), x, m)), y=0, n+1), y, n);
end:

##Tie sequences: seq(HowManyGames(A,i,i),i=0..10)
##Soccer, A={1}
#1, 2, 6, 20, 70, 252, 924, 3432, 12870, 48620, 184756
#A000984, Central binomial coefficients
##Old-time Basketball, A={1,2}
#1, 2, 14, 84, 556, 3736, 25612, 177688, 1244398, 8777612, 62271384
#A036692, sequence refers to Zeilberger and basketball
##Basketball, A={1,2,3}
#1, 2, 14, 106, 784, 6040, 47332, 375196, 3001966, 24190148, 196034522
#A122680, sequence refers to Zeilberger and basketball
##American football, A={2,3,6,7,8}
#1, 0, 2, 2, 6, 24, 62, 206, 712, 2076, 7284
#Not in OEIS