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

Help:=proc():
print(`UmbT(P,var,N),`):
end:

#1. Write a procedure, UmbT(P,var,N), that inputs a polynomial P in a list of variables,
# var, say [x1,x2,x3, ..xk], and outputs the first N terms of the sequence that is the
# "umbral transform" under the umbra P^n under the umbra
# x1^a1 .... xk^ak -> (a1+...+ak)!/(a1!*...*ak!)
UmbT:=proc(P,var,N) local n:
[seq(Umb(P^n,var),n=1..N)]:
end:

(* Outputs:

a) UmbT(x + y, [x, y], 20);
[2, 6, 20, 70, 252, 924, 3432, 12870, 48620, 184756, 705432, 
  2704156, 10400600, 40116600, 155117520, 601080390, 2333606220, 
  9075135300, 35345263800, 137846528820]

A000984		Central binomial coefficients: binomial(2*n,n) = (2*n)!/(n!)^2.

b) UmbT(x+y+z,[x,y,z],20);
[3, 15, 93, 639, 4653, 35169, 272835, 2157759, 17319837, 140668065, 1153462995, 9533639025,
 79326566595, 663835030335, 5582724468093, 47152425626559, 399769750195965, 3400775573443089,
 29016970072920387, 248256043372999089]

A002893		a(n) = Sum_{k=0..n} binomial(n,k)^2 * binomial(2*k,k).

c) UmbT(x1 + x2 + x3 + x4, [x1, x2, x3, x4], 20);
[4, 28, 256, 2716, 31504, 387136, 4951552, 65218204, 878536624, 
  12046924528, 167595457792, 2359613230144, 33557651538688, 
  481365424895488, 6956365106016256, 101181938814289564, 
  1480129751586116848, 21761706991570726096, 
  321401321741959062016, 4766118425002290943216]

A002895		Domb numbers: number of 2n-step polygons on diamond lattice.

d) UmbT(1 + 2*x + 3*y, [x, y], 20);
[6, 48, 432, 4104, 40176, 400896, 4053888, 41396832, 425922624, 
  4408432128, 45850157568, 478784058624, 5016552777216, 
  52713659768832, 555293062213632, 5862253175834112, 
  62006462404125696, 656969858722013184, 6971258566157672448, 
  74074268834613817344]

Couldn't find this one in OEIS

e) UmbT((1 + x)*(1 + y)*(1 + z), [x, y, z], 20);
[16, 550, 24136, 1176106, 60777016, 3261771136, 179738342800, 
  10099540746586, 576074600766736, 33252430325188900, 
  1938057963369265600, 113864747946771636400, 
  6735066708925272520000, 400678521748749771325600, 
  23955811161103318460282656, 1438502844846334412424107866, 
  86710074762372112298341225600, 5244447235864729051389952560556, 
  318158612150147855918832514776496, 
  19353875636742012377756474564054476]

Couldn't find this one, either.

*)

#2. Have a fix: The issue is how IsCond3G deals with ties, which is not an issue for the odd terms.
# Made the following modification:

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

(* Now NuC behaves as expected, returning half the actual number.

> Zs(10);
     [0, 0, 12, 0, 540, 180, 21000, 13440, 785820, 711900]

> 2*[seq(NuC(v, [[0, 1, 1], [0, 0, 1]]), v = 1 .. 10)];
     [0, 0, 12, 0, 540, 180, 21000, 13440, 785820, 711900]


*)

#3. Done.

#### Included without change from C18.txt ####
#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:





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