#Homework 23 by Ben Miles for Experimental math OK to post #1 PT:=proc() local S,L,S1: L:=[3]: S:=[]: while nops(S)<92 do S:=[op(S),SplitUp(L)]: L:=JC1(L): od: S: end: #2 Mat:=proc() local a,i,j,table,co,E1,E2,E3: table:=PT(): for i from 1 to 92 do for j from 1 to 92 do co:=0: E1:=SplitUp(table[i]): E2:=table[j]: for E3 in E1 do if E3= E2 then co:=co+1: fi: od: a[i,j]:=co: od: od: Matrix(a): end: #4 RLT1:=proc(S,N) local P,u,s: P:=1: for s in S do P:=P+u^s: od: print(P): LIF(P,u,N): end: #5 RLT2:=proc(S,N) local P,u,s,m,P1: P1:=0: for s in S do P1:=P1+u^s: od: m:=max(S): P:=add(u^i,i=0..m)-P1: print(P): LIF(P,u,N): end: ##################### #LIF: 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: #Counting unlabelled rooted trees #A[n]=number of rooted unalabelled trees with n vertices #A[1]=1, A[2]=1, A[3]=2, #Removing the root, leads to a MULTISET of smaller rooted trees #HOW MANY with one vertex? 1/(1-x)^A[1] #HOW MANY with two vertices 1/(1-x^2)^A[2] #... #How MANY with i vertices? 1/(1-x^i)^A[i] #A[1]*x+A[2]*x^2+A[3]*x^3+A[4]*x^4+...= #x/(1-x)^A[1]/(1-x^2)^A[2]/(1-x^3)^A[3]/..... #URT(N): the first N terms of the enumerating #sequence for the numnber of UNLABELLED rooted trees #with n vertices URT:=proc(N) local A,f,n,x,i: A:=[1]: for n from 2 to N do f:=x*mul(1/(1-x^i)^A[i],i=1..nops(A)): f:=taylor(f,x=0,n+1): A:=[op(A),coeff(f,x,n)]: od: A: 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: SplitUp:=proc(L) local a,L1,L2,i: a:=nops(L): for i from 1 to a-1 do: L1:=L[1..i]: L2:=L[i+1..a]: if split(L1,L2) then RETURN([op(SplitUp(L1)),op(SplitUp(L2))]): fi: od: [L]: end: Help:=proc(): print(`JC(w), YN(INI,n), JC1(w) `): 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: