# John Miller # hw23.txt - Homework assigned Thu, Apr 12, due Apr 15. # It is OK to post this. Help := proc(): print(`PT(), Mat(), EstConwayConst()`): print(`RLT1(S,N), RLT2(S,N)`): print(`SplitUp(L), split(L,R), LargestAB(L)`): print(`JC(w), YN(INI,n), JC1(w), YN1(INI,n)`): print(`LIF(P,u,N)`): end: ###################################################################### # Question #1 # PT(): Outputs the set of Conway's 92 elements. PT := proc() local PeriodicTable, Compound::list, i, NewAtoms, AtomsToCheck, OldAtom: Compound := [3]: PeriodicTable := {[3]}: AtomsToCheck := {[3]}: while AtomsToCheck <> {} do: OldAtom := AtomsToCheck[1]: Compound := JC1(OldAtom): NewAtoms := convert(SplitUp(Compound),set) minus PeriodicTable: PeriodicTable := PeriodicTable union NewAtoms: AtomsToCheck := AtomsToCheck union NewAtoms minus {OldAtom}: od: PeriodicTable: end: # end of PT() ###################################################################### # Question #2 # Mat(): Outputs the 0-1 matrix such that the entry (i,j) has a 1 iff # the "molecule" JC1(w), applied to atom i, has atom j in it. Mat := proc() local PTable::list, i, j, A::array, sizePT, splitMolecule: PTable := convert(PT(),list): sizePT := nops(PTable): A := Array(1..sizePT,1..sizePT): for i from 1 to sizePT do: splitMolecule := SplitUp(JC1(PTable[i])): for j from 1 to sizePT do: if evalb(PTable[j] in splitMolecule) then A[i,j] := 1: else A[i,j] := 0: fi: od: od: A: end: # end of Mat(): ###################################################################### # Question 3 # EstConwayConst(): Use the LinearAlgebra package to find the largest eignvalue # of Mat(). EstConwayConst := proc() local M::Matrix, EigValues, x: with(LinearAlgebra): Digits := 20: M := convert(Mat(),Matrix): #EigValues := evalf(Eigenvalues(M,output='list')): #LargestAB(EigValues): CharacteristicPolynomial(M,x): end: # end of EstConwayConst() # The output I get is 1.30299... which is close to Conway's constant # 1.30357... but still different enough that I wonder why I don't get # a better match. ###################################################################### # Question 4 # RLT1(S,N): 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. # The generating function u is given by: # u = t*(1 + add(u^a/a!, a in S)) = t*P RLT1 := proc(S::set, N::posint) local a, P, u, n, egf_coeff: P := 1 + add(u^a/a!, a in S): egf_coeff := LIF(P,u,N): [seq(egf_coeff[n]*n!,n=1..N)]: end: # end of RLT1(S,N) # Some example outputs: # RLT1({1},20) outputs: # 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800, 39916800, 479001600, 6227020800, 87178291200, 1307674368000, 20922789888000, 355687428096000, 6402373705728000, 121645100408832000, 2432902008176640000 # which is OEIS A000142 (n!) # RLT1({1,2},20) outputs: # 1, 2, 9, 60, 540, 6120, 83790, 1345680, 24811920, 516650400, 11992503600, 307069963200, 8598348158400, 261387760233600, 8573572885878000, 301809119163552000, 11349727401396384000, 454104511068656448000, 19261139319649202976000, 863322072620761353600000 # which is OEIS A036774 (# of labeled rooted binary trees) # RLT1({1,2,3},20) outputs: # 1, 2, 9, 64, 620, 7620, 113610, 1992480, 40194000, 916927200, 23341071600, 655922836800, 20169411662400, 673645440468000, 24285190867938000, 939899116892736000, 38870133445791648000, 1710655202853140544000, 79826043011286892320000, 3936948118406837614080000 # which is OEIS A036775 # RLT1({1,3},20) outputs: # 1, 2, 6, 28, 200, 1920, 22260, 299040, 4596480, 80035200, 1558972800, 33556723200, 790404014400, 20220183984000, 558377964144000, 16556553965952000, 524646641046528000, 17693631589810176000, 632750850059328000000, 23917010504518010880000 # is not in OEIS. ###################################################################### # Question 4 # RLT2(S,N): 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 # belong to S. # The generating function u is given by: # u = t*(1 + add(u^a/a!, a not in S)) # = t*(exp(u) - add(u^a/a!, a in S)) = t*P RLT2 := proc(S::set, N::posint) local a, P, u, n, egf_coeff: P := exp(u) - add(u^a/a!, a in S): egf_coeff := LIF(P,u,N): [seq(egf_coeff[n]*n!,n=1..N)]: end: # end of RLT2(S,N) # Some example outputs: # RLT2({1},20) outputs: # 1, 0, 3, 4, 65, 306, 4207, 38424, 573057, 7753510, 134046671, 2353898196, 47602871329, 1013794852266, 23751106404495, 590663769125296, 15806094859299329, 448284980183376078, 13515502344669830287, 430050399901131972780 # which is OEIS A060356 # RLT2({1,2},20) outputs: # 1, 0, 0, 4, 5, 6, 427, 1968, 6561, 220510, 2129171, 13847736, 337904437, 5156062926, 54298310445, 1192150218496, 24147409593089, 364887230459454, 8145781717395223, 197451127561855320 # which is not in OEIS. # RLT2({1,2,3},20) outputs: # 1, 0, 0, 0, 5, 6, 7, 8, 2529, 11350, 36971, 104556, 10182757, 99054970, 632882265, 3303250096, 165364954369, 2689602118254, 28186612549255, 233809699635780 # which is not in OEIS. # RLT({1,3},20) outptus: # 1, 0, 3, 0, 65, 6, 3787, 1184, 427905, 286750, 79563671, 92894616, 22050773329, 39814849270, 8526143835315, 22030869781696, 4386817826347009, 15360562101209598, 2898552816857784799, 13208608054252514040 # which is not in OEIS. ###################################################################### #LIF(P,u,N): Inputs an expression P of the variable #u that possesses #a Maclaurin expansion in u, and outputs the #first N terms of the series expansion of the #unique f.p.s u(t) satisfying the functional eq. #u(t)=t*P(u(t)). For example, #LIF(1+u^2,u,10); should return something with #Catalan numbers LIF:=proc(P,u,N) local n: [seq(coeff(taylor(P^n, u=0,n),u,n-1)/n,n=1..N)]: end: ###################################################################### # LargestAB(L): Inputs a list of complex numbers, and find the one with # largest absolute value. # For example, LargestAB([5,-3,-4,6,-7,1]); should give -7. LargestAB := proc(L::list) local i, absVal, maxAB, maxIndex: maxAB := 0: maxIndex := 1: for i from 1 to nops(L) do absVal := abs(evalf(L[i])): if absVal > maxAB then maxAB := absVal: maxIndex := i: fi: od: L[maxIndex]: end: # end of LargestAB(L) ###################################################################### # SplitUp(L): Inputs a list L in {1,2,3} and outputs it as a list # of lists [L1,L2,L3,..], where L = [op(L1),op(L2),op(L3),...] and # L1, L2, L3 are atoms. SplitUp := proc(L::list) option remember: local i, word1::list, word2::list: if nops(L) = 1 then return([L]): fi: for i from 1 to nops(L)-1 do: word1 := L[1..i]: word2 := L[i+1..nops(L)]: if split(word1,word2) then return([word1,op(SplitUp(word2))]): fi: od: [L]: end: # end of SplitUp(L) ############################# stuff from C22.txt ################################### #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: ##################################################################333 #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: ####################################################################3 #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: #########################################################################3 #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: ###################################################################### 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: