#OK to post homework
#Karnaa Mistry, 10/11/20, Assignment HW #9

with(combinat):


# 1. 

# (i)
# DiagSeq2(1/(1-1*x-9*y+1*x*y-2*x^3+1*y^3),x,y,50)[51] = 656121696907513857036622746341010698644640606651318753632076909371800565836900539

# (ii)
# DiagWalks3D({[1,2,2],[9,0,0],[0,2,2]},30)[31] = 433160




# 2.

#WordsMod(a,S,x): inputs a positive integer a, subsest S of {1,2,..., a} and a variable x and outputs the rational function that is the
#generating function of finite sequence of positive integers that leave remainder that belongs to S when divided by a
#Since S cannot have an element 0, we'll say that for some integer p, performing (p*a mod a) would yield "a", instead of zero.
WordsMod:=proc(a,S,x) local f,g,h,i:

f:=0:

for i from 1 to a do
    if (i in S) then
        f:=f+x^i:
    fi:
od:

g:=(1-x^a):

h:=f/g:

RETURN(factor(normal(1/(1 - h)))):

end:

# seq(coeff(taylor(WordsMod(5,{1,2,3,4},x),x=0,101),x,i),i=0..100)[101] = 118272656870560296827634782280 ways




# 3.

#ResComps(S,x): inputs a finite set S of positive integers, and outputs the generating function for a(n) defined as the number
#of compositions of n where each component can be any positive integer except members of S
ResComps:=proc(S,x) local i,f,s,L,T:

f:=x/(1-x):

f:= f - add(x^s,s in S):

RETURN(factor(normal(1/(1-f)))):

end:


# seq(coeff(taylor(ResComps({2,3},x),x=0,301),x,i),i=0..300)[301] = 179789662603291780480882832586094549526183607154633616127794901




# 4.

# (i) is in the OEIS, A00984.  f(n) = (4n+2)/n * f(n-1), f(0) = 1
# (ii) is in the OEIS, A001850.  f(n) = (6n-3)/n*f(n-1) + (1-n)/n*f(n-2), f(0) = 1, f(1) = 3
# (iii) is in the OEIS, A036692.
#    {(352*n^3 + 2768*n^2 + 7088*n + 5920)*f(n + 1) + (-484*n^3 - 4048*n^2 - 11184*n - 10204)*f(n + 2) + (-352*n^3 - 3120*n^2 - 9114*n - 8766)*f(n + 3)
#    + (55*n^3 + 515*n^2 + 1570*n + 1560)*f(n + 4), f(0) = 1, f(1) = 2, f(2) = 14, f(3) = 84} (base listtorec from Maple, unsure of simplification)
# (iv) is in the OEIS, A192365.  The Maple gfun of this fails.




# 5. 

# (i) is in the OEIS, A006480.  f(n) = (27n^2+27n+6)/(n^2+2n+1) * f(n-1), f(0) = 1
# (ii) is in the OEIS, A268542.
#    {(-567*n^3 - 4104*n^2 - 9549*n - 7140)*f(n + 1) + (-1281*n^3 - 10553*n^2 - 28546*n - 25400)*f(n + 2) + (-672*n^3 - 6208*n^2 - 18752*n
#     - 18400)*f(n + 3) + (42*n^3 + 430*n^2 + 1424*n + 1504)*f(n + 4), f(0) = 1, f(1) = 4, f(2) = 42, f(3) = 520} (base listtorec from Maple)
# (iii) is in the OEIS, A112019.
#    {(-59*n^3 - 448*n^2 - 1084*n - 848)*f(n + 1) + (-295*n^3 - 2535*n^2 - 7185*n - 6740)*f(n + 2) + (-2301*n^3 - 22074*n^2 - 69941*n - 73044)*f(n + 3) 
#    + (118*n^3 + 1250*n^2 + 4336*n + 4896)*f(n + 4), f(0) = 1, f(1) = 5, f(2) = 55, f(3) = 749} (base listtorec from Maple)




# 6. 

#DiagSeq4(f,x1,x2,x3,x4,N): 4-dimensional analog of DiagSeq3
DiagSeq4:=proc(f,x1,x2,x3,x4,N) local i:

if subs({x1=0,x2=0,x3=0,x4=0},denom(f))=0 then
 print(`The denominator of`, f, `should have a non-zero constant term `):
 RETURN(FAIL):
fi:

[seq(coeff(taylor(coeff(taylor(coeff(taylor(coeff(taylor(f,x1=0,N+1),x1,i),x2=0,N+1),x2,i),x3=0,N+1),x3,i),x4=0,N+1),x4,i),i=0..N)]:

end:


#DiagWalks4D(S,N): 4-dimensional analog of DiagWalks3D
DiagWalks4D:=proc(S,N) local s,x,y,z,w:
if not (type(S, set) and type(N,integer) and N>=0) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

DiagSeq4(1/(1-add(x^s[1]*y^s[2]*z^s[3]*w^s[4], s in S)),x,y,z,w,N ) :

end: