#OK to post homework
#Blair Seidler, 2021-03-28 Assignment 16

with(combinat):

Help:=proc():
print(`ComputCons(F, P, n, d)`):
end:

# 1.

#ComputCons(F, P, n, d): that inputs an F, that is a product and/or quotient of expressions
#of the form (a*n+b)! (e.g. binomial(2*n,n)/10^n), a polynomial P in n, a variable name n,
#and a positive integer d, and outputs the value of
#Sum(((-1)^n*P(n)*F(n),n=0..infinity) to d decimal digits.
ComputCons:=proc(F, P, n, d) local total,fcur,pcur,fratio,i,tol:
Digits:=d:
fratio:=simplify(expand(subs(n=n+1,F)/F)):
total:=1:
tol:=10^(-d-1):
fcur:=evalf(subs(n=0,F)):
pcur:=evalf(subs(n=0,P)):
total:=fcur*pcur:
for i from 1 while fcur*pcur>tol do
  fcur:=evalf(fcur*subs(n=i-1,fratio)):
  pcur:=evalf(subs(n=i,P)):
  total:=total+(-1)^(i mod 2)*fcur*pcur:
od:
evalf(total):
end:

#### Analysis ####
# I ran this procedure and the naive infinite sum for k=1..100 and d=10000
# ComputCons ran in 319.5s, but the Sum took 2826.9s

# 2. 
# I was unable to run LT(7,7) using the code from class. I tried running it overnight and it
# didn't finish.
# The values of LT(i,i) for i=1..6 are 1,0,1,0,4,236, which does not appear in the OEIS.

# 3.
# I ran seq(nops(LT(i,2)),i=2..8)
# 0, 1, 4, 19, 112, 771, 6088
#A127548   O.g.f.: Sum_{n>=0} n!*(x/(1+x)^2)^n.

# 4.

#I've spent some time thinking about this one, but I don't have a proof yet. I think that
#it should be related to derangements somehow, since what we need is a second row which is
#a permutation of n with no conflicts. The only difference from derangements is that instead
#of trying to avoid pi(i)=i for 1..n, we are trying to avoid pi(i) in {i,i+1} for i=1..n-1.
#We are really ignoring pi(n) because it is the unused number in the second row.

#### From C16.txt ####

#IsGood1(T): Given a potential trapezoid (in the right-angled format)
# that you know that its horiz. ,  vertical 
#and its NW->SE diagonals for the truncated version (where you chopped
#the last row), returns true iff the new one is also OK
#Try:
#IsGood1([[1,2,3],[2,1]]);
IsGood1:=proc(T) local n,k,i:
k:=nops(T):
n:=nops(T[1]):
for i from 1 to nops(T[k]) do
#checking that all NW->SW diagonals are OK
 
 if member(T[k][i],{seq(T[k-j][i+j],j=1..k-1)}) then
   RETURN(false):
 fi:

#Checking that all vertical lines are distinct
if member(T[k][i],{seq(T[j][i],j=1..k-1)}) then
   RETURN(false):
 fi:
od:

true:

end:

#LT(n,k): The SET OF REDUCED n by k LATIN TRAPEZOIDS
#The bottom line is 1,...,n, and each floor above it
#of length n-1,n-2, .., n-k+1 has all DIFFERENT entries
#also VERTICALLY all different also A[i][i],A[i+1][i+1], ... are
#all different
LT:=proc(n,k) local S, Sold,PNF,old1,f,new1:
option remember:

if k=1 then
  RETURN({[[seq(i,i=1..n)]]}):
fi:

Sold:=LT(n,k-1):


PNF:=permute(n,n-k+1):
S:={}:
for old1 in Sold do
  for f in PNF do
    new1:=[op(old1),f]:


   if IsGood1(new1) then
      S:=S union {new1}:
   fi:
  od:
od:
S:
end: