> #ok to post
;
> #Yifan Zhang, 10/8/2020, hw9
;
> 
;
> with(combinat);
[Chi, bell, binomial, cartprod, character, choose, composition, conjpart, 
decodepart, encodepart, eulerian1, eulerian2, fibonacci, firstcomb, firstpart,
firstperm, graycode, inttovec, lastcomb, lastpart, lastperm, multinomial, 
nextcomb, nextpart, nextperm, numbcomb, numbcomp, numbpart, numbperm, partition
, permute, powerset, prevcomb, prevpart, prevperm, randcomb, randpart, randperm
, rankcomb, rankperm, setpartition, stirling1, stirling2, subsets, unrankcomb,
unrankperm, vectoint]

;
> #Q1.
;
> #1. DiagSeq2(1/(1-1*x-7*y+1*x*y-1*x^3+9*y^3),x,y,49)[50]
;
> #333077213928786617477722190187481245846499973315721325230123564513744107
;
> #2. {[1,9,9], [7,1,1], [1,9,9]}
;
> #DiagWalks3D({[1,9,9], [7,1,1], [1,9,9]},30);
> #[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
;
> 
;
> 
;
> #Q2.
;
> WordsMod:= proc(a,S,x) local s:
> option remember:
> 
> # input S is a subset of {1,2,3,4,....,a}
> #WordsMod(5,{1,2,3,4}, x)
> 
> factor( 1/(1-add(x^s, s in S)/(1-x^a))):
> end:
> 
;
> GFseq:=proc(f,x,N) local  f1,i:
> f1:=taylor(f,x=0,N+1):
> [seq(coeff(f1,x,i),i=0..N)]:
> end:
> WordsMod(5,{1,2,3,4,5},x);
(x-1)/(2*x-1)

;
> GFseq((x-1)/(2*x-1),x,99)[100];
316912650057057350374175801344

;
> #316912650057057350374175801344
;
> 
;
> 
;
> 
;
> #Q3.
;
> #Input a finite set S of positive integers, output the generating function for a(n): nunber of composition of n where each component cannot be in S.
;
> #Such as RecComps({2,3},x)
;
> #
;
> ResComps:= proc(S,x) local s:
> option remember:
> 
> if not (type(S,set) and min(op(S))>0 ) then
>  print(`Bad input`):
>  RETURN(FAIL):
> fi:
> 
> factor(1/(1-add(x^s,s in S ))):

> 
> end:
> ResComps({2,3},x)
-1/(x^3+x^2-1)

;
> GFseq(-1/(x^3+x^2-1),x,299)[300];
1346693590052167920449786597761583494

;
> 
;
> 
;
> 
;
> #Q4.
;
> #
;
> #
;
> #
;
> #DiagWalks2D(S,N): Inputs  a set of "atomic steps" in the 2-dimenional square lattice, with positve steps, of the form [a,b]
> #For example for a forward moving Knight the set if {[2,1],[1,2]} for a forward moving King it  is {[1,0],[0,1],[1,1]}
> #outputs the first N+1 terms of the sequence
> #Number of walks from [0,0] to [n,n] using the atomic steps. For example try:
> #DigSeq({[1,2],[2,1]},40);
> DiagWalks2D:=proc(S,N) local s,x,y:
> if not (type(S, set) and type(N,integer) and N>=0) then
>  print(`Bad input`):
>  RETURN(FAIL):
> fi:
> 
> DiagSeq2(1/(1-add(x^s[1]*y^s[2], s in S)),x,y,N  ) :
> end:
> #DiagSeq2(f,x,y,N): Given a rational function f of the variables x and y such that the denominator does not vanish at (0,0)
> #outputs the list whose i-th entry is the coefficient of x^(i-1)*y^(i-1) in the 2-variable power-series expansion of f
> #Try:
> #DiagSeq2(1/(1-x-y),x,y);
> DiagSeq2:=proc(f,x,y,N) local i:
> 
> if subs({x=0,y=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(f,x=0,N+1),x,i),y=0,N+1),y,i),i=0..N)]:
> 
> end:
> #1. A984
;
> DiagWalks2D({[0,1],[1,0]},40)
[1, 2, 6, 20, 70, 252, 924, 3432, 12870, 48620, 184756, 705432, 2704156, 
10400600, 40116600, 155117520, 601080390, 2333606220, 9075135300, 35345263800,
137846528820, 538257874440, 2104098963720, 8233430727600, 32247603683100, 
126410606437752, 495918532948104, 1946939425648112, 7648690600760440, 
30067266499541040, 118264581564861424, 465428353255261088, 1832624140942590534,
7219428434016265740, 28453041475240576740, 112186277816662845432, 
442512540276836779204, 1746130564335626209832, 6892620648693261354600, 
27217014869199032015600, 107507208733336176461620]

;
> gfun[listtorec](%,f(n))
[{(-n-2)*f(n+2)+(6+4*n)*f(n+1), f(0) = 1, f(1) = 2}, ogf]

;
> 
;
> #Guess: f(n) =((6+4*n)/(n+2))*f(n-1), f(0) = 1, f(1) = 2
;
> 
;
> #2. 
;
> DiagWalks2D({[0,1],[1,0],[1,1]},40)
[1, 3, 13, 63, 321, 1683, 8989, 48639, 265729, 1462563, 8097453, 45046719, 
251595969, 1409933619, 7923848253, 44642381823, 252055236609, 1425834724419, 
8079317057869, 45849429914943, 260543813797441, 1482376214227923, 
8443414161166173, 48141245001931263, 274738209148561921, 1569245074591690083, 
8970232353223635949, 51313576749006450879, 293733710358893793729, 
1682471873186160624243, 9642641465118083682429, 55294491352291112747007, 
317241780630136241094657, 1820991621200098527441027, 10457362855894862001750669
, 60078868458555230983696959, 345299757825442889707393857, 
1985346154034162284608274707, 11419126147845924397833900957, 
65701922725618214591910684159, 378150244155138145169182750209]

;
> gfun[listtorec](%,f(n))
[{(-n-3)*f(n+3)+(15+6*n)*f(n+2)+(-n-2)*f(n+1), f(0) = 1, f(1) = 3, f(2) = 13},
ogf]

;
> #guess: f(n) = ((15+6*n)*f(n-1)-(n+2)*f(n-2))/(n+3), f(0) = 1, f(1) = 3, f(2) = 13
;
> 
;
> 
;
> #3. A36682
;
> DiagWalks2D({[0,1],[1,0],[0,2],[2,0]},5)
[1, 2, 14, 84, 556, 3736]

;
> gfun[listtorec](%,f(n))
[{(29*n+116)*f(n+4)+(-418-120*n)*f(n+3)+(-2002-773*n)*f(n+2)-678*f(n+1), f(0) =
1, f(1) = 2, f(2) = 14, f(3) = 84}, ogf]

;
> #guess: f(n) = ((352*n^3+2768*n^2+7088*n+5920)*f(n-3)+(-484*n^3-4048*n^2-11184*n-10204)*f(n-2)+(-352*n^3-3120*n^2-9114*n-8766)*f(n-1))/(-1*55*n^3+515*n^2+1570*n+1560)), f(0) = 1, f(1) = 2, f(2) = 14, f(3) = 84
;
> 
;
> 
;
> #4
;
> DiagWalks2D({[0,1],[0,2],[1,0],[2,0],[1,1],[2,2]},40)
[1, 3, 22, 165, 1327, 10950, 92045, 783579, 6733966, 58294401, 507579829, 
4440544722, 39000863629, 343677908223, 3037104558574, 26904952725061, 
238854984979423, 2124492829796598, 18927927904130617, 168888613467092895, 
1508973226894216106, 13498652154574126523, 120886709687492946083, 
1083687170568092836350, 9723660300694365146989, 87322023363899721467343, 
784797155029928933462614, 7058384549327774062817985, 63525027490219905080597647
, 572078727948766828614679830, 5154896021503193637126695629, 
46475144470335949808346622587, 419221380265186873217492981134, 
3783331574173309950816296187025, 34158702327707968075053323650645, 
308541107632605100591602320772210, 2788040608584603973862476702488265, 
25202873888846971990310549726981859, 227906767043212151853904258820045110, 
2061638450059846741700344736532466833, 18655567955631232569396589476504942877]

;
> gfun[listtorec](%,f(n))
FAIL

;
> 
;
> 
;
> 
;
> #Q5
;
> #1.
;
> DiagSeq3:=proc(f,x,y,z,N) local i:
> 
> if subs({x=0,y=0,z=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(f,x=0,N+1),x,i),y=0,N+1),y,i),z=0,N+1),z,i),i=0..N)]:
> 
> end:
> DiagWalks3D:=proc(S,N) local s,x,y,z:
> if not (type(S, set) and type(N,integer) and N>=0) then
>  print(`Bad input`):
>  RETURN(FAIL):
> fi:
> 
> DiagSeq3(1/(1-add(x^s[1]*y^s[2]*z^s[3], s in S)),x,y,z,N  ) :
> end:
> DiagWalks3D({[0,0,1],[0,1,0],[0,0,1]},50)
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

;
> DiagWalks3D({[0,0,1],[0,1,0],[0,0,1],[1,1,0],[1,0,1],[0,1,1]},50)
[1, 4, 42, 520, 7090, 102144, 1525776, 23380368, 365130810, 5786380600, 
92774019052, 1501646797248, 24498046138384, 402329384914240, 6645072333486720,
110293868867458080, 1838511122725436250, 30762545845461663240, 
516457722792721132500, 8696582711157975572400, 146835895664138204448540, 
2485266624268853914959360, 42157397527828247726954880, 
716556351190046311597161600, 12201942968211237198350281200, 
208134928903175159702319457344, 3555830679926229005075061257376, 
60836819329493358733292740546048, 1042259677840638259492432048065280, 
17878405805956405560805535960396160, 307035478496101645580861028622654336, 
5278643189711854927868638533509575744, 90844825001410805657489064172778794842,
1564930724096710135146829959035675638440, 
26982512386206869455023952107027146825700, 
465628368604516908867820871730799473867408, 
8041657743459510320462629847980612873329892, 
138989457278403173476481541185549849571771136, 
2403979582946502330471619673933968184168916960, 
41607939250187986849488309498865757934964797280, 
720614911745512461611918130596030003676289336300, 
12488149132570126719418439010195646125397100011920, 
216544860562054593134160906994219210611674235418920, 
3756996029385637316936199664140175184307538677619200, 
65217708880091867222809412770976203671887863352995200, 
1132692240424335329461008559322860625660096843564044800, 
19682044848299254204780187424606161624293922138158848000, 
342161666360328010168145138304650840941581198111307680000, 
5950957372613238042937365323128808212313545207994095090800, 
103544899011058176270363371934019286365853202810957930168000, 
1802392915688776002461574283464036323570366183012829246320352]

;
> gfun[listtorec](%,f(n))
[{(-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}, ogf]

;
> #Guess: f(n) =( -(-672*n^3-6208*n^2-18752*n-18400)*f(n-1)-(-1281*n^3-10553*n^2-28546*n-25400)*f(n-2)-(-567*n^3-4104*n^2-9549*n-7140)*f(n-3))/(42*n^3+430*n^2+1424*n+1504), f(0) = 1, f(1) = 4, f(2) = 42, f(3) = 520
;
> 
;
> DiagWalks3D({[0,0,1],[0,1,0],[0,0,1],[1,1,0],[1,0,1],[0,1,1],[1,1,1]},50);
[1, 5, 55, 749, 11251, 178835, 2949115, 49906925, 860905315, 15071939255, 
266982872905, 4774722189275, 86070844191775, 1561948324845095, 
28507384046515555, 522867506128197869, 9631571375362268515, 
178094411589895650815, 3304192479145474141741, 61487420580006795749999, 
1147310054439368264752621, 21460298038569915209597885, 
402303281700710751346492045, 7557003072139605981578257499, 
142216452439487217855645948751, 2680953577121013488247925591235, 
50618646171303180157424116235125, 957111527948218432660657321864175, 
18121739540270261728632231230077015, 343543557236065454833612395025716335, 
6520361879189996271307536840076027315, 123890073917796829080870582463811238125,
2356385395407636535251297627384060933475, 
44861516384576265010060080056738042597135, 
854859275252481703527187912584359902505365, 
16303683385055931306937559135507039289051575, 
311190647408479503515035135447889843225653825, 
5944262823858514467535608131557811648002762625, 
113627280554332441572175257060408538344973808425, 
2173522796137552287361594832873182554514084408175, 
41603303651437038812509444164890710254979329256765, 
796818898586340848384595439852318646653791885982025, 
15270271969281638238724672178615586077051065833360775, 
292804205106533694230148888112886916503795367402916085, 
5617461573225955169248923401253208251816762860322509065, 
107826322191169954469150116083770101402658959404009097505, 
2070721317569227785240110965397887711335116126944031387125, 
39785206672880588731577259303187378878045846931120125328475, 
764744172323001796507565473412652088019110691985844385750095, 
14706094191942557902358875596276810126406906869637116414501755, 
282916070222353435645433592830184350982846989381549672223125905]

;
> gfun[listtorec](%,f(n))
[{(-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}, ogf]

;
> #Guess: f(n) = (-(-2301*n^3-22074*n^2-69941*n-73044)*f(n-1)-(-295*n^3-2535*n^2-7185*n-6740)*f(n-2)-(-59*n^3-448*n^2-1084*n-848)*f(n-3))/(118*n^3+1250*n^2+4336*n+4896), f(0) = 1, f(1) = 5, f(2) = 55, f(3) = 749
;
> 
;
> 
;
> 
;
> 
;
> #Q6.
;
> DiagSeq4:=proc(f,x,y,z,k,N) local i:
> 
> if subs({x=0,y=0,z=0,k=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,x=0,N+1),x,i),y=0,N+1),y,i),z=0,N+1),z,i),k=0, N+1),k,i), i=0..N)]:
> 
> end:
> DiagWalks4D:=proc(S,N) local s,x,y,z,k:
> 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]*k[4], s in S)),x,y,z,k,N  ) :
> end:
> 
;
> 
;
>  
;