#!/usr/local/bin/maple
#-*- maplev -*-

# Nathaniel Shar
# HW 22
# Experimental Mathematics

# It is okay to link to this assignment on the course webpage.

Help := proc(): print(`JC1(w), PT(), DetectUnknownAtoms(L, elements), Mat(), ConwayConstant(), LIF(P,u,N), RLT1(S, N), RLT2(S, N)`): end:

# From C22:

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:

#############
# Problem 1 #
#############

PT := proc() local elements, todo, s, products: option remember:
    todo := {[3]}:
    elements := {}:
    while todo <> {} do:
        s := todo[1]:
        todo := todo minus {s}:
        products := DetectUnknownAtoms(JC1(s), elements):
        todo := todo union products:
        elements := elements union products:
    od:
    return elements:
end:

DetectUnknownAtoms := proc(L, elements) local atoms, a, unknown:
    unknown := {}:
    atoms := SplitUp(L):
    for a in atoms do:
        if not member(a, elements) then:
            unknown := unknown union {a}:
        fi:
    od:
    return unknown:
end:

#############
# Problem 2 #
#############

Mat := proc() local elts: option remember:
    elts := PT():
    return Matrix(92, 92, ((i,j) -> numboccur(elts[j], SplitUp(JC1(elts[i]))))):
end:


#############
# Problem 3 #
#############

with(LinearAlgebra):

ConwayConstant := proc() local L:
    L := Eigenvalues(Mat(), output = `list`):
    return map(attributes, sort([seq](setattribute(abs(evalf(ev)), evalf(ev)), ev=L), `>`))[1]:
end:

# returns 1.30357726903430

#############
# Problem 4 #
#############

# From c23:

LIF := proc(P,u,N) local n:
    [seq(coeff(taylor(P^n, u=0,n),u,n-1)/n,n=1..N)]:
end:

RLT1 := proc(S, N) local i, egf:
    egf := LIF(1 + add(u^j/j!, j=S), u, N):
    return seq(i!*egf[i], i=1..N):
end:

# In OEIS: {1}, {1,2}, {1,2,3}
# Not in OEIS: {1,3}

#############
# Problem 5 #
#############

RLT2 := proc(S, N) local i, egf:
    egf := LIF(exp(u) - add(u^j/j!, j=S), u, N):
    return seq(i!*egf[i], i=1..N):
end:

# In OEIS: {1} (combinatorial interpretation not mentioned, though)
# Not in OEIS: {1,2}, {1,2,3}, {1,3}