# This will compute the dimension, using Ngai and Wang's method.
# This will assume that the IFS has maps $f_i(x) = 1/beta*x + b_i, 
# where the DS = {b_1, b_2, ..., b_n}.

with(linalg):
with(LinearAlgebra);
ComputeDim := proc(DS, beta)
    
    local kappa, d1, d2, CL, NL, NT, A, NS, ns, DI, i, n, NS2, NS3, 
          NTO, NI, c, NT2, j, NID, x, P, Poly, NT2C;
    # CL - CurrentLevel
    # NL - NextLevel
    # NT - NeighbourhoodTree
    # A  - AdjacenyGraph
    # NS - NeighbourhoodSet
    # DI - DigitIndex
    # NI - NextItemToLookAt
    # NID - NextItemToLookAtDone;
 
    for i from 1 to nops(DS) do 
       DI[DS[i]] := i;
    od;

    # First thing I want to do is compute kappa.
    kappa := 0;
    for d1 in DS do for d2 in DS do
        if evalf(abs(d1-d2)>kappa) then kappa := abs(d1-d2) fi;
    od; od;

    kappa := evalf(kappa / (beta-1));

    NI := {{{0}}};
    NT := {};
    NID := {};

    while NI <> {} do 
 
        CL := NI[1];
        NID := NID union {CL};
        CL := CL[1];
        NI := NI[2..-1];

        lprint(nops(NI), nops(NID));
        NL := {};
        for d1 in DS do 
            for c in CL do 
                NL := {[c*beta+d1, c, DI[d1]]} union NL;
            od;
        od;
        NL := expand(NL);
        NL := simplify(NL);

        NS := {};
        for d1 in DS do
            ns := select(y->y[1]=d1, NL);
            ns := [op(ns)];
            ns := sort(ns, (y,z)->y[-1]<z[-1]);
            NS := NS union {ns[1]};
        od;
 
        NL := map(y->y[1], NL);
        NS := map(y->[GetNeighbours(y[1], NL, kappa),
                      GetNeighbours(y[2], CL, kappa),
                      y[3]], NS);
        NS := map(y->[{y[1]}, {y[2]}, y[3]], NS);
        
        NI := NI union (map(y->y[1], NS) minus NID);
        NT := NT union NS;
    od;
  
    NT2 := map(y->op(y[1..2]), {op(NT)});
    for i from 1 to nops(NT2) do 
        NT2C[NT2[i]] := i;
    od;

    lprint("Making a matrix, size "||(nops(NT2))|| " by "
                                  ||(nops(NT2)));
    A := Matrix(nops(NT2), nops(NT2), 0);
    for i in NT do 
        A[NT2C[i[1]],NT2C[i[2]]] := A[NT2C[i[1]],NT2C[i[2]]] +1;
    od;

    Poly := MinimalPolynomial(A, x);
    P := max(fsolve(Poly));
    Poly := factor(Poly);
    if not type(Poly, `+`) then 
        Poly := [op(Poly)];
        for p in Poly do 
            if abs(subs(x=P, p)) < 10^(-Digits/2) then 
                Poly := p;
                break;
            fi;
        od;
    fi;

    RETURN(P, Poly, evalf(log(P)/log(beta)), Log(P)/Log(beta));
end:
    
GetNeighbours := proc(n, NL, kappa)
    local ns;
    ns := select(y->evalf(Re(evalf(abs(y-n)))<evalf(kappa)), NL);
    ns := map(y->y-n, ns);
    RETURN(ns);
end:

Digits := 20;

# 3-gaskets of order infinity.
K := ComputeDim({seq(exp(2*Pi*I*i/3),i=0..2)}, 2, 2);

# 3-gasket of order 2
tau := RootOf(x^2-x-1, index=1);
K := ComputeDim({seq(exp(2*Pi*I*i/3),i=0..2)}, tau, 3);

j := 3:
N := 6:
T[N] := [seq(exp(2*Pi*I*i/6),i=1..N)]:
K := ceil(N/4)-1:
tf[N] := (2+2*add(cos(2*Pi*i/N),i=1..K)):
t[N] := RootOf((tf[N]-1)-tf[N]*la^(j-1)+la^j,index=2):

# 6-Gasket of order infinity
ComputeDim(T[6], tf[6], 6);

# 6-Gasket of order 2
ComputeDim(T[6], t[6], 6);

j := 3:
N := 8:
T[N] := [seq(exp(2*Pi*I*i/N),i=1..N)]:
K := ceil(N/4)-1:
tf[N] := (2+2*add(cos(2*Pi*i/N),i=1..K)):
t[N] := RootOf((tf[N]-1)-tf[N]*la^(j-1)+la^j,index=2):

# 8-Gaskets of order infinity
TT := ComputeDim(T[8], tf[8]); factor(TT[2]);

# 8-Gaskets of order 2
TT := ComputeDim(T[8], t[8]); factor(TT[2]);

N := 12:
T[N] := [seq(exp(2*Pi*I*i/N),i=1..N)]:
K := ceil(N/4)-1:
tf[N] := (2+2*add(cos(2*Pi*i/N),i=1..K)):
j := 2;
t[N, j] := RootOf((tf[N]-1)-tf[N]*la^(j-1)+la^j,index=2):
j := 3;
t[N, j] := RootOf((tf[N]-1)-tf[N]*la^(j-1)+la^j,index=2):

# 12-Gasket of order infinity
TT := ComputeDim(T[12], tf[12]); factor(TT[2]);
# 12-Gasket of order 1
TT := ComputeDim(T[12], t[12, 2]); factor(TT[2]);
# 12-Gasket of order 2
TT := ComputeDim(T[12], t[12, 3]); factor(TT[2]);