# 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: