# Please do not post homework
# Lucy Martinez, 04-03-22, Assignment 18
read `C18.txt`:
with(combinat):
#----------------------Problem 2-----------------------#
# Using procedure NuC(v,G) with G=[[0,1,1],[0,0,1]] write a one-line procedure SeqCoStupid(N)
# that inputs a positive integer N and outputs the first N terms of the sequence "Number of voting-profiles 
# (NOT vote-count profiles) with 2*n-1 voters that lead to 1->2->3->1"

SeqCoStupid:= proc(N) local n:
[seq(NuC(2*n-1,[[0,1,1],[0,0,1]]),n=1..N)]:
end:
SeqCoStupid(7);
#        [0, 6, 270, 10500, 392910, 14478156, 529404876]


#----------------------Problem 3-----------------------#
# [A little bit challenging, do your best]---Done in class with Dr. Z.
# Read this gem and implement Theorem 1 to write a procedure SeqCoSmart(N) that does the same thing as
# SeqCoStupid(N), but hopefully faster.
# [Hint, use the taylor command to exapand the rational function in t to N+1 terms, extract the first N 
# coefficients of t, getting polynomials in the 6 variables, x123, ..., x321, and then apply the 
# "umbra" as described in the paper. How much faster is SeqCoSmart(20) than SeqCoStupid(20)?

# Umb1(M,var): inputs a monomial M in the list of variables var and applies the "umbra" 
# c*var[1]^a1*...*var[k]^ak -> c*(a1+...+ak)!/(a1!*...*ak!)
Umb1:=proc(M,var) local f,d,i,c:
 for i from 1 to nops(var) do
  d[i]:=degree(M,var[i]):
 od:
c:=normal(M/mul(var[i]^d[i],i=1..nops(var))):
 
if {seq(normal(diff(c,var[i])),i=1..nops(var))}<>{0} then
 RETURN(FAIL):
fi:
 
c*add(d[i],i=1..nops(var))!/mul(d[i]!,i=1..nops(var));
end:
 
# Umb(P,var): applying the above umbra to the polynomial P
Umb:= proc(P,var) local P1,i:
if P=0 then
 RETURN(0):
fi:
P1:=expand(P):
add(Umb1(op(i,P1),var),i=1..nops(P1)):
end:
 
SeqCoSmart:=proc(N) local f, var,f1,i,L :
f:=(1+t^3*x123*x231*x312)*(t^3*x312*x123*x231)/((1-t^2*x123*x321)*(1-t^2*x213*x312)*(1-t^2*x312*x123)*(1-t^2*x132*x231)*(1-t^2*x123*x231)*(1-t^2*x231*x312)):

f1:=taylor(f,t,N+4):
var:=[x123,x132,x213,x231,x312,x321]:

L:=[seq(coeff(f1,t,i),i=1..N)]:

[seq(Umb(L[i],var),i=1..nops(L))]:
end:
2*SeqCoSmart(7);
#              [0, 0, 6, 0, 540, 6, 21000]

SeqCoStupid(7);
#        [0, 6, 270, 10500, 392910, 14478156, 529404876]

time(SeqCoStupid(20));
#                        933.718
time(SeqCoSmart(20));
#                        0.093
# Since time(SeqCoStupid(20)), time(SeqCoSmart(20))=0.093, then SeqCoSmart(20) is 933.625 times faster


#----------------------Problem 4-----------------------#
# Try to estimate the limiting probability of the Condorcet scenario, by cranking out as many terms as you can
# and dividing by 6^(2*n-1) (and at the end mulitplying by 2)


#----------------------Problem 5-----------------------#
# [Optional challenge, I don't know the answer, 10 dollars, and possibly a paper]
# Find a nice combinatorial proof (if possible a bijective one) of the fact that we conjectured in class
# (and that can be made into an ugly proof) that the number of vote-count-Conrorcet scenariors with 2*v-1
# voters is binomial(v+3,5) (recall that the total number of vote-count scenarios is binomial(2*v+4,5)


#----------------------Problem 6-----------------------#
# A sequence of integers is a quasi-polynomial of degree d and period p if for i=1, ..,p the subsequence 
# {L[p*n+i]}, is a polynomial of degree d in n. Using GP(L,n), as a subroutine, first write a procedure
# GQP1(L,p,n) that inputs a list L and tries to guess a quasi-polynomial of period p that fits the data, 
# and then 
# write GQP(L,n) that tries out p=1, p=2, until either a success or failure. The output should be a list of 
# polynomials in n,let's call it P, such that if n mod p=i (but we take i=1,...p, rather than the usual
# i=0,..,p-1) then L[n]=P[i][n]. For example GQP([1,-2,3,-4,5,-6,7,-8,9,-10,11,-12,13,-14]) should output [n,-n]

GQP1:=proc(L,p,n) local i,j,m,Li:
for i from 1 to p do
 Li:=[seq(L[p*m+i],m=0..floor((nops(L)-i)/p))]: 
 # print(Li): this would print out each subsequence
 print(GP(Li,n)):
od:
end:

GQP1([1,-2,3,-4,5,-6,7,-8,9,-10,11,-12,13,-14,15,-16],4,n);
#                    4 n - 3
#                   -4 n + 2
#                    4 n - 1
#                   -4 n