> #Weiji Zheng, 11/01/2020
;
> #Assignment 16
;
> 
;
> 
;
> #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-2)*p_3(n-3))
;
> #       = p_3(n-1) + (n-1)*p_3(n-2) + (n^2-3*n+2)*p_3(n-3)
;
> #     0 = 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)
;
> 
;
> #ope(n,N)=N^3 - N^2 - (n + 2)*N - (n + 1)(n + 2)
;
> #order of recurrence ==> 3
;
> 
;
> #code.
;
> #SeqFromRec(ope,n,N,Ini,K): Given the first L-1
> #terms of the sequence Ini=[f(1), ..., f(L-1)]
> #satisfied by the recurrence ope(n,N)f(n)=0
> #extends it to the first K values
> 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)
ope := N^3-N^2-(n+2)*N-(n+2)*(n+1)
;
> 
;
> SeqFromRec(ope,n,N,[1,1,1],200)[200] 
499792773985469608207915734551449152758014145434000190057722878158099799204935545287271844300939648649245434108606684301694686331807780102161764048092100159116900670693897587287578850120870208045036789702743874406143710190279405832776711541973801348099997696
;
> 
;
> 
;
> #Q3
;
> #Write a procedure 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 Inputs a positive integer n;Outputs the sum of binomial(n,k)^6 from k=0 to n
;
> 
;
> S6 := proc(n) local i:
> add(binomial(n,i)^6,i=0..n):
> end:
> 
;
> #S6seq Inputs a positive integer N; Outputs the first N terms of the sequence S6(n)
;
> 
;
> S6seq:=proc(N) local n:
> seq(S6(n),n=1..N):
> end:
> 
;
> #Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
> #poly coffs.
> #of maximum DEGREE+ORDER<=MaxC
> #e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
> Findrec:=proc(f,n,N,MaxC)
> local DEGREE, ORDER,ope,L:
> 
> for L from 0 to MaxC do
> for ORDER from 0 to L do
>  DEGREE:=L-ORDER:
> if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
>  print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
>  RETURN(FAIL):
> fi:
>  ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):
>      if ope<>FAIL then
>        RETURN(ope):
>       fi:
>  od:
> od:
> FAIL:
> 
> end:
> 

;