#OK to post homework
#Blair Seidler, 4/3/22, Assignment 18
with(combinat):
with(ListTools):

Help:=proc():
print(`SeqCoStupid(N), SeqCoSmart(N)`):
end:

#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:

#3. SeqCoSmart(N): 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.
#
# Stole the Umbra code from C19.txt and modified Zs from C19.txt

SeqCoSmart:=proc(N) local f,var,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));
f:=taylor(f,t,2*N+2):
var:=[x123,x132,x213,x231,x312,x321]:

L:=[seq(expand(coeff(f,t,2*i-1)),i=1..N)]:

[seq(Umb(L[i],var),i=1..nops(L))]:
end:

(* How much faster is SeqCoSmart(20) than SeqCoStupid(20)?

A lot faster:

time(SeqCoStupid(20));
                            750.531
time(SeqCoSmart(20));
                             11.718

*)

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

(* [0., 0.1111111111, 0.1388888889, 0.1500342936, 0.1559523129, 0.1596281531, 0.1621368284,
    0.1639591603, 0.1653432359, 0.1664302878, 0.1673067051, 0.1680283199, 0.1686328275,
    0.1691466001, 0.1695886429, 0.1699729990, 0.1703102698, 0.1706086088, 0.1708743895,
    0.1711126661, 0.1713274981, 0.1715221835, 0.1716994294, 0.1718614788, 0.1720102057,
    0.1721471882, 0.1722737646, 0.1723910772, 0.1725001073, 0.1726017020, 0.1726965969,
    0.1727854335, 0.1728687741, 0.1729471134, 0.1730208885, 0.1730904870, 0.1731562538,
    0.1732184968, 0.1732774919, 0.1733334868, 0.1733867048, 0.1734373473, 0.1734855969,
    0.1735316194, 0.1735755654, 0.1736175724, 0.1736577659, 0.1736962609, 0.1737331628, 0.1737685685]

Best guess for what this is approching (from Inverse Symbolic Calculator is:
sum(1/(binomial(2*n, n)/(1/2*n^3 + 7/2*n^2 + 10*n - 11)), n = 1 .. infinity)
but I haven't really thought about why this might be.
*)

#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]
#GP1(L,n,d): Inputs a list of numbers L and a variable n and a non-neg. d outputs the polynomial of degree d
#P(n), such that P(i)=L[i], for all i from 1 to nops(L). OR FAIL
#nops(L)>=d+2
GP1i:=proc(L,n,d,IL) local a,i,P,eq,var:
if nops(IL)<d+2 then
 ERROR(`List too short`):
 fi:
#P=a[0]+a[1]*n+...+a[d]*n^d
P:=add(a[IL[i+1]]*n^i,i=0..d):
#print(P):
var:={seq(a[IL[i+1]],i=0..d)}:
eq:={seq(L[IL[i]]=subs(n=IL[i],P),i=1..d+2)}:

#print(eq,var):

var:=solve(eq,var):
#print(var):
 if var=NULL then
   RETURN(FAIL):
 fi:
P:=subs(var,P):
if {seq(L[IL[i]]-subs(n=IL[i],P),i=d+3..nops(IL))}<>{0} then
  RETURN(FAIL):
 fi:
P:
end:

#GPi(L,n,IL): Inputs a list of numbers L and a variable n and outputs the polynomial of degree <=nops(L)-2
#P(n), such that P(i)=L[i], for all i from 1 to nops(L). OR FAIL
#nops(L)>=d+2
GPi:=proc(L,n,IL) local d,P:
for d from 0 to nops(IL)-2 do
#print(L,n,d,IL):
 P:=GP1i(L,n,d,IL):
 if P<>FAIL then
    RETURN(P):
 fi:
od:
FAIL:
end:

GQP1:=proc(L,p,n) local i,j,IL,PL,P:
PL:=[]:
for i from 1 to p do
  IL:=[seq(p*j+i,j=0..floor((nops(L)-i)/p))]:
#print(L,n,IL):
  P:=GPi(L,n,IL):
  if P=FAIL then RETURN(FAIL): fi:
  PL:=[op(PL),P]:
od:
PL:
end:


GQP:=proc(L,n) local p,P:
for p from 1 to floor(nops(L)/4) do
 P:=GQP1(L,p,n):
 if P<>FAIL then
    RETURN(P):
 fi:
od:
FAIL:
end:

(* Notes:

I couldn't figure out how to do this with the original GP, because the n values are off when
I used a subarray. So I modified GP and GP1 to take an additional argument, IL, which is the
list of positions to consider within L.

A couple of sample outputs:

> GQP([1, -2, 3, -4, 5, -6, 7, -8, 9, -10, 11, -12, 13, -14], n);

                            [n, -n]



> L := [seq(cos(Pi*i/2)*i^2 + cos(Pi . i)*i^3, i = 1 .. 40)];

 L := [-1, 4, -27, 80, -125, 180, -343, 576, -729, 900, -1331, 
        1872, -2197, 2548, -3375, 4352, -4913, 5508, -6859, 8400, 
       -9261, 10164, -12167, 14400, -15625, 16900, -19683, 22736, 
      -24389, 26100, -29791, 33792, -35937, 38148, -42875, 47952, 
      -50653, 53428, -59319, 65600]

> GQP(L, n);
                  [-n^3 , n^3  - n^2 , -n^3 , n^3  + n^2 ]

*)


#### Included from C19.txt ####
#C19.txt: March 31, 2021. Maple code for Lecture 19, Using the umbral approach to count Condorcet vote-profile scenarios
Help19:=proc(): print(`Umb1(M,var), Umb(P,var), Zs(N) `): end:
#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
#Debugged after class
Umb:=proc(P,var) local P1,i:
if P=0 then
 RETURN(0):
fi:

P1:=expand(P):
if not type(P1,`+`) then
RETURN(Umb1(P,var)):
fi:


add(Umb1(op(i,P1),var),i=1..nops(P1)):
end:


#### Included from C18.txt ####

#C18.txt: Maple code for Lecture 18, original and Generaized Condorcet scenario
Help18:=proc(): print(`IsCond3G(C,G), Cond3G(v,G) , SeqC(N,G), SeqCe(N,G), `):
print(`SeqCo(N,G), NuC(v,G) `):end:

#i beats j if the number of voters that prefer i to j is larger than
#the number of voters that prefer j to i
#
#S[i,j]:=Number of voters that prefer i to j 
#i beats j if S[i,j]-S[j,i]>=1, 1 is the cut-off
#
#Let G be a 2 by 3 cut-off matrix, [[0,G12,G13],[0,0,G23]]

#IsCond3G(C,G): inputs a vote-count profile outputs true iff 1->2->3->1,
IsCond3G:=proc(C,G) local i,P,T:
P:=permute(3):

for i from 1 to nops(C) do
 T[P[i]]:=C[i]:
od:

#1->2->3->1
if (T[[1,2,3]]+T[[1,3,2]]+T[[3,1,2]])- (T[[2,1,3]]+T[[2,3,1]]+T[[3,2,1]])>=G[1][2] and
(T[[1,2,3]]+T[[2,1,3]]+T[[2,3,1]])-(T[[1,3,2]]+T[[3,1,2]]+T[[3,2,1]])>=G[2][3] and
-(T[[2,3,1]]+T[[3,2,1]]+T[[3,1,2]])+(T[[2,1,3]]+T[[1,2,3]]+T[[1,3,2]])<G[1][3] then
 true:
 else
 false:
fi:
end:


Cond3G:=proc(v,G) local A,a,C:
A:=Comps(v,6):
C:={}:
#C is the SET of vote-count profiles with v voters that lead to 1->2->3->1
for a in A do
if IsCond3G(a,G) then
 C:=C union {a}:
fi:
od:
C:
end:

#Nuc(v,G): The total number of voter-profiles (not to be confused with vote-count-profiles) with 1->2->3->1 with v voters
NuC:=proc(v,G) local A,su,a,i:
A:=Comps(v,6):
su:=0:
for a in A do
 if IsCond3G(a,G) then
  su:=su+add(a)!/mul(a[i]!,i=1..nops(a)):
 fi:
od:
su:
end:





#SeqCo(N,G): The first N terms of the sequence
#Number of G-Condorcet vote-COUNT-profiles with 2*n-1 voters
SeqCo:=proc(N,G) local n:
[ seq(nops(Cond3G(2*n-1,G)),n=1..N)]:
end:

#SeqCe(N,G): The first N terms of the sequence
#Number of G-Condorcet vote-COUNT-profiles with 2*n voters
SeqCe:=proc(N,G) local n:
[ seq(nops(Cond3G(2*n,G)),n=1..N)]:
end:

#SeqC(N,G): The first N terms of the sequence
#Number of G-Condorcet vote-COUNT-profiles with v voters
SeqC:=proc(N,G) local n:
[ seq(nops(Cond3G(n,G)),n=1..N)]:
end:

Help17:=proc(): print(`Comps(n,k), GP1(L,n,d), GP(L,n) , IsCond3(C) `):end:



#Comps(n,k): The set of all compositions of n into exactly k parts (with non-negative entrties)
Comps:=proc(n,k) local S,a,S1,s:
option remember:
if k=0 then
 if n=0 then
   RETURN({[]}):
 else
  RETURN({}):
 fi:
fi:
S:={}:

for a from 0 to n do
 S1:=Comps(n-a,k-1):
 S:=S union {seq([op(s),a],s in S1)}:
 od:
 S:
end:



#GP1(L,n,d): Inputs a list of numbers L and a variable n and a non-neg. d outputs the polynomial of degree d
#P(n), such that P(i)=L[i], for all i from 1 to nops(L). OR FAIL
#nops(L)>=d+2
GP1:=proc(L,n,d) local a,i,P,eq,var:
if nops(L)<d+2 then
 ERROR(`List too short`):
 fi:
#P=a[0]+a[1]*n+...+a[d]*n^d
P:=add(a[i]*n^i,i=0..d):
var:={seq(a[i],i=0..d)}:
eq:={seq(L[i]=subs(n=i,P),i=1..d+2)}:
var:=solve(eq,var):
 if var=NULL then
   RETURN(FAIL):
 fi:
P:=subs(var,P):
if {seq(L[i]-subs(n=i,P),i=d+3..nops(L))}<>{0} then
  RETURN(FAIL):
 fi:
P:
end:




#GP(L,n): Inputs a list of numbers L and a variable n and outputs the polynomial of degree <=nops(L)-2
#P(n), such that P(i)=L[i], for all i from 1 to nops(L). OR FAIL
#nops(L)>=d+2
GP:=proc(L,n) local d,P:
for d from 0 to nops(L)-2 do
 P:=GP1(L,n,d):
 if P<>FAIL then
    RETURN(P):
 fi:
od:
FAIL:
end:


#IsCond3(C): inputs a vote-count profile outputs true iff 1->2->3->1,
IsCond3:=proc(C) local i,v,P,T:
P:=permute(3):
v:=add(C):
for i from 1 to nops(C) do
 T[P[i]]:=C[i]:
od:

#1->2->3->1
if T[[1,2,3]]+T[[1,3,2]]+T[[3,1,2]]>v/2 and 
T[[1,2,3]]+T[[2,1,3]]+T[[2,3,1]]>v/2 and
T[[2,3,1]]+T[[3,2,1]]+T[[3,1,2]]>v/2 then
 true:
 else
 false:
fi:
end: