# Alex Gentul Homework 23 OK to post info1:=proc() print(`PT(), Mat()`); end: info1(); # Write a procedure PT() (PT for Periodic Table), that inputs nothing, and outputs, the set of Conway's 92 elements. Verify that # indeed you got 92 of them. # Hint: Start with [3], and keep applying SplitUp(L), of the previous homework assignment, that is based on the true/false # procedure # split(w1,w2) # that you should use as a black box. PT:=proc() local L,L1,i; L:=[3]; L1:=[]; while nops(L)<92 do L:=[op(L),SplitUp(L1)]; L1:=JC1(L); end; L; end: # Write a procedure Mat(), that inputs nothing, and outputs the 0-1 matrix such that entry (i,j) has a 1 if and only if the # "molecule" (that may be an atom) # JC1(w) # applied to atom i, has atom j in it. Mat:=proc() local M,i,j,A1,A2,a,n,m; M:=PT(); for i from 1 to 92 do for j from 1 to 92 do n:=0; A1:=SplitUp(M[i]); A2:=M[j]; for a in A2 do if a = A2 then n:=n+1; end; end; m[i,j]=n; end; Matrix(m); end: # Using LIF(P,u,N), and exponential generating functionolgy, write a procedure, # RLT1(S,N), # that inputs a set of positive integers S, and a positive integer N, and outputs the list of size N whose i-th entry is the # number of rooted labelled trees on i vertices where each vertex either has outdegree 0 (a leaf) or an outdegree that belongs to # S. RLT1:=proc(S,N) local P,u,n; P:=1; for n in S do P:=P+u^n; end; LIF(P,u,N); end: # Using LIF(P,u,N), and exponential generating functionolgy, write a procedure, # RLT2(S,N) , # that inputs a set of positive integers S, and a positive integer N, and outputs the list of size N whose i-th entry is the # number of rooted labelled trees on i vertices where each vertex either has outdegree 0 (a leaf) or an outdegree that DOES NOT # belongs S. RLT2:=proc(S,N) local P1,P,u,n,m,i; P1:=0; for n in S do P1:=P1+u^n; end; m:=max(S); P:=add(u^i,i=0..m)-P1; LIF(P,u,N); end: ########## other stuff ########### SplitUp:=proc(L) local L1,L2,i; for i from 1 to nops(L)-1 do L1:=L[1..i]; L2:=L[i+1..nops(L)]; if split(L1,L2) then return [op(SplitUp(L1)),op(SplitUp(L2))]; end; end; [L]; end: split:=proc(L,R) local n,m: if nops(L)=0 or nops(R)=0 then RETURN(1): fi: n:=op(nops(L),L): m:=op(1,R): if n=m then RETURN(false): fi: if n>=4 and m<=3 then RETURN(true): fi: if n=2 then if m=1 and nops(R)>2 and op(2,R)<>1 and op(3,R)<>op(2,R) then RETURN(true): fi: if m=1 and nops(R)=2 and op(2,R)<>1 then RETURN(true): fi: if m=1 and nops(R)>=3 and op(2,R)=1 and op(3,R)=1 then RETURN(true): fi: if m=3 and nops(R)>=2 and op(2,R)<>3 and not(nops(R)>=4 and op(2,R)=op(3,R) and op(3,R)=op(4,R)) then RETURN(true): fi: if m=3 and nops(R)=1 then RETURN(true): fi: if m>=4 then RETURN(true): fi: fi: if n<>2 then if nops(R)>=5 and op(1,R)=2 and op(2,R)=2 and op(3,R)=1 and op(4,R)<>1 and op(5,R)<>op(4,R) then RETURN(true): fi: if nops(R)=4 and op(1,R)=2 and op(2,R)=2 and op(3,R)=1 and op(4,R)<>1 then RETURN(true): fi: if nops(R)>=5 and op(1,R)=2 and op(2,R)=2 and op(3,R)=1 and op(4,R)=1 and op(5,R)=1 then RETURN(true): fi: if nops(true)>=4 and op(1,R)=2 and op(2,R)=2 and op(3,R)=3 and op(4,R)<>op(3,R) and not(nops(R)>=6 and op(4,R)=op(5,R) and op(5,R)=op(6,R)) then RETURN(true): fi: if nops(R)=2 and op(1,R)=2 and op(2,R)=2 then RETURN(true): fi: if nops(R)=3 and op(1,R)=2 and op(2,R)=2 and op(3,R)>=4 then RETURN(true): fi: if nops(R)>=3 and op(1,R)=2 and op(2,R)=2 and op(3,R)>=4 and op(4,R)<>op(3,R) then RETURN(true): fi: fi: false: end: #The Conway Look-And-Say sequence #JC(w): inputs any list in and outputs its #Yusra transform #a1^c1,a2^c2,..., ar^cr ->[a1,c1,a2,c2,a3,c3, ...] JC:=proc(w) local i,a1,c1,w1,J1: option remember: if w=[] then RETURN([]): fi: a1:=w[1]: for i from 1 to nops(w) while w[i]=a1 do od: c1:=i-1: w1:=w[c1+1..nops(w)]: J1:=JC(w1): [c1,a1,op(J1)]: end: #YN(INI,n): the first n terms of the #sequence : "length of JC^n(INI) YN:=proc(INI,n) local s,i,Cu: s:=[]: Cu:=INI: for i from 1 to n do s:=[op(s), nops(Cu)]: Cu:=JC(Cu): od: s: end: #JC1(w): inputs any list in and outputs its #Yusra transform, not using recursion #a1^c1,a2^c2,..., ar^cr ->[a1,c1,a2,c2,a3,c3, ...] JC1:=proc(w) local i,a1,c1,w1,J1,L: option remember: if w=[] then RETURN([]): fi: L:=[]: w1:=w: while w1<>[] do a1:=w1[1]: for i from 1 to nops(w1) while w1[i]=a1 do od: c1:=i-1: L:=[op(L),c1,a1]: w1:=w1[c1+1..nops(w1)]: od: L: end: #YN1(INI,n): the first n terms of the #sequence : "length of JC^n(INI) YN1:=proc(INI,n) local s,i,Cu: s:=[]: Cu:=INI: for i from 1 to n do s:=[op(s), nops(Cu)]: Cu:=JC1(Cu): od: s: end: LIF:=proc(P,u,N) local n: [seq(coeff(taylor(P^n, u=0,n),u,n-1)/n,n=1..N)]: end: #GuessH1(L,ORD,DEG,n,N): Inputs a sequence of numbers L #and pos. integer ORD and non-neg. integer DEG and letters #n and N and outputs an linear recurrence operator with #poly. coeffs. (conj.) annihilating the sequence L #nops(L)>=(1+ORD)*(1+DEG)+ORD+6 GuessH1:=proc(L,ORD,DEG,n,N) local ope,a,i,j,var,eq,n1,var1: if nops(L)<(1+ORD)*(1+DEG)+ORD+6 then print(`We need at least`, (1+ORD)*(1+DEG)+ORD+6, `terms in L`): RETURN(FAIL): fi: ope:=add(add(a[i,j]*n^i*N^j,i=0..DEG),j=0..ORD): var:={seq(seq(a[i,j],i=0..DEG),j=0..ORD)}: eq:={seq(add(subs(n=n1, coeff(ope,N,j))*L[n1+j],j=0..ORD), n1=1..(1+DEG)*(1+ORD)+6)}: var1:=solve(eq,var): ope:=subs(var1,ope): if ope=0 then RETURN(FAIL): fi: ope:=subs({seq(v=1, v in var)}, ope): ope:=add(factor(coeff(ope,N,j))*N^j, j=0..degree(ope,N)): end: