with(combinat);

pet_cycleind_symm :=
proc(n)
local p, s;
option remember;

    if n=0 then return 1; fi;

    expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;

pet_flatten_term :=
proc(varp)
local terml, d, cf, v;

    terml := [];

    cf := varp;
    for v in indets(varp) do
        d := degree(varp, v);
        terml := [op(terml), seq(v, k=1..d)];
        cf := cf/v^d;
    od;

    [cf, terml];
end;

cycles_prod :=
proc(cyca, cycb)
local ca, cb, lena, lenb, res, vlcm;

    res := 1;

    for ca in cyca do
        lena := op(1, ca);

        for cb in cycb do
            lenb := op(1, cb);

            vlcm := lcm(lena, lenb);
            res := res*a[vlcm]^(lena*lenb/vlcm);
        od;
    od;

    res;
end;


automaton :=
proc(N, M, K)
option remember;
local idx_slots, idx_cols, idx_syms, res, a, b, c, sim, flat_sim,
    sym, flat_sym, flat_a, flat_b, flat_c,
    cyc_a, cyc_b, len_a, len_b, p, q;

    if N > 1 then
        idx_slots := pet_cycleind_symm(N);
    else
        idx_slots := [a[1]];
    fi;

    if M > 1 then
        idx_cols := pet_cycleind_symm(M);
    else
        idx_cols := [a[1]];
    fi;

    if K > 1 then
        idx_syms := pet_cycleind_symm(K);
    else
        idx_syms := [a[1]];
    fi;


    res := 0;

    for a in idx_slots do
        flat_a := pet_flatten_term(a);
        for b in idx_cols do
            flat_b := pet_flatten_term(b);

            sim := cycles_prod(flat_a[2], flat_b[2]);
            flat_sim := pet_flatten_term(sim);

            for c in idx_syms do
                flat_c := pet_flatten_term(c);

                sym := cycles_prod(flat_b[2], flat_c[2]);
                flat_sym := pet_flatten_term(sym);

                p := 1;

                for cyc_a in flat_sim[2] do
                    len_a := op(1, cyc_a);
                    q := 0;

                    for cyc_b in flat_sym[2] do
                        len_b := op(1, cyc_b);

                        if len_a mod len_b = 0 then
                            q := q + len_b;
                        fi;
                    od;

                    p := p*q;
                od;

                res := res +
                p*flat_a[1]*flat_b[1]*flat_c[1];
            od;
        od;
    od;

    res;
end;