> #ok to post 
;
> #Yifan Zhang, hw16, 10/28/2020
;
> 
;
> #Q1.
;
> #An involution is a permutation whose cycle structure consists only of cycles of length 1 and 2. Let p_3(n) be the number of permutations whose cycle #structure consists of cycles of length 1, 2, or 3. Find a linear recurrence satisfied by p_3(n), and express the linear recurrence operator ope(n,N) annihilating  it. What is the order of the recurrence? Using procedure SeqFromRec in in ComboProject4.txt find p_3(200).
> #p_3(n)=p_3(n-1)+(n-1)*p_3(n-2)-(n-1)*(n-2)*p_3(n-3)
;
> #p_3(n)-p_3(n-1)-(n-1)*p_3(n-2)-(n-1)*(n-2)*p_3(n-3)=0
;
> #p_3(n+3)-p_3(n+2)-(n+2)*p_3(n+1)-(n+2)*(n+1)*p_3(n)=0
> #-(n+2)*(n+1)*p_3(n)-(n+2)*p_3(n+1)-p_3(n+2)+p_3(n+3)=0
;
> #Linear Recurrence. 
;
> SeqFromRec:=proc(ope,n,N,Ini,K)
> local ope1,gu,L,n1,j1:
> ope1:=Yafe(ope,N)[2]:
> L:=degree(ope1,N):
> if nops(Ini)<>L then
>  ERROR(`Ini should be of length`, L):
> fi:
> 
> ope1:=expand(subs(n=n-L,ope1)/N^L):
> 
> gu:=Ini:
> 
> for n1 from nops(Ini)+1 to K do
> gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
> j1=1..L)]:
> od:
> 
> gu:
> 
> end:
> Yafe:=proc(ope,N) local i,ope1,coe1,L: 
> if ope=0 then
>  RETURN(1,0):
> fi:
> ope1:=expand(ope):
> L:=degree(ope1,N):
> coe1:=coeff(ope1,N,L):
> ope1:=normal(ope1/coe1):
> ope1:=normal(ope1):
> ope1:=
> convert(
> [seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
> factor(coe1),ope1:
> end:
> ope:=N^3-N^2-(n+2)*N-(n+2)*(n+1)
Typesetting:-mprintslash([(ope := N^3-N^2-(n+2)*N-(n+2)*(n+1))],[N^3-N^2-(n+2)*
N-(n+2)*(n+1)])

;
> SeqFromRec(ope,n,N,[1,1,1],10) 
[1, 1, 1, 10, 26, 96, 552, 2316, 12108, 72696]

;
> SeqFromRec(ope,n,N,[1,1,1],200)[200] 
4997927739854696082079157345514491527580141454340001900577228781580997992049355\
45287271844

;
> #Q2.
;
> #[Optional Challenge] Let p_k(n) be the number of permutations whose cycle structure consists of cycles of length ≤ k. Find a linear recurrence satisfied by p_k(n), and express the linear recurrence operator ope(n,N) annihilating it. What is the order of the recurrence?
;
> #for length of 1,2, or 3
#ope:=N^3-N^2-(n+2)*N-(n+2)*(n+1)
#for length of 1,2,3 or 4
#ope:=N^4-N^3-(n+3)*N^2-(n+3)*(n+2)*N-(n+3)*(n+2)*(n+1)
#for length of 1,2,3,4, or 5
#ope:=N^5-N^4-(n+4)*N^3-(n+4)*(n+3)*N^2-(n+4)*(n+3)*(n+2)*N-(n+4)*(n+3)*(n+2)*(n+1)
;
> #ope(n,N) := N^k-N^(k-1)-(n+k-1)*N^(k-2)-(n+k-1)*(n+k-2)*N^(k-3)-...-(n+k-1)*(n+k-2)*...*(n+k-(k-2))*N-(n+k-1)*(n+k-2)*...*(n+2)*(n+1)
;
> #the order of recurrence is k
;
> 
;
> #Q3.
;
> #Write a procedure Sum4(n) that inputs a positive integer n and outputs the sum of binomial(n,k)^6 from k=0 to n. Call is S6(n). Then write a one-line procedure S6seq(N) that inputs a positive integer N and outputs the first N terms of the sequence S6(n).
> #Using procedure Findrec in ComboProject4.txt find a linear recurrence operator annihilating it (equivalently a recurrence) and the initial condition, and write procedure S6seqClever(N) that uses procedure SeqFromRec from the same Maple package to do the same thing as S6seq(N).
> 
> #What are
> 
> #time(S6seq(1000));
> 
> #and
> 
> #time(S6seqClever(1000)); ?
> 

;
> S6:=proc(n) local k, S:
> S:=0:
> for k from 0 to n do
> S:= S+binomial(n,k)^6;
> od:
> S:
> end:
> S6(3)
1460

;
> S6seq:=proc(N):
> [seq(S6(i), i=1..N)];
> end:
> S6seq(10)
[2, 66, 1460, 54850, 2031252, 86874564, 3848298792, 180295263810, 8709958973540
, 433617084579316]

;
> 
;
> time(S6seq(1000))
8.079

;
> S6seqClever:=proc(m) local ope, n, N:
> ope:=-8*(6*n+5)*(2*n+3)*(510578*n+1212419)*(6*n+7)*(n+1)^3/(n+3)/(68821*n+158376)/(n
> +4)^5+1/3*(1220155462*n^7+16171613287*n^6+95508822435*n^5+325758743385*n^4+
> 689826913083*n^3+900310603788*n^2+665180159400*n+213076106100)/(n+3)/(68821*n+
> 158376)/(n+4)^5*N-1/3*(329726969*n^7+6818567734*n^6+59948741865*n^5+
> 290983523330*n^4+843114141824*n^3+1459404674232*n^2+1398114248708*n+
> 572049603720)/(n+3)/(68821*n+158376)/(n+4)^5*N^2-2/3*(4178086*n^6+81090768*n^5+
> 657966606*n^4+2852370092*n^3+6957071718*n^2+9037771995*n+4877455590)/(68821*n+
> 158376)/(n+4)^5*N^3+N^4:
> SeqFromRec(ope, n, N, [2, 66, 1460, 54850], m):
> 
> end:
> S6seqClever(10)
[2, 66, 1460, 54850, 2031252, 86874564, 3848298792, 180295263810, 8709958973540
, 433617084579316]

;
> time(S6seqClever(1000))
.221

;
> #Conclusion: S6seqClever() is much much faster than S6seq().
;
> 
;
> 
;
> 
;
> 
;
> 
;
> 
;
> 
;