pet_cycleind_symm :=
proc(n)
option remember;
    local l;

    if n=0 then return 1; fi;

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

pet_cycleind_edg :=
proc(n)
option remember;
local all, term, termvars, res, l1, l2, inst1, u, v,
    uidx, vidx;

    if n=0 or n=1 then return 1; fi;

    all := 0:
    for term in pet_cycleind_symm(n) do
        termvars := indets(term); res := 1;

        # edges on different cycles of different sizes
        for uidx to nops(termvars) do
            u := op(uidx, termvars);
            l1 := op(1, u);

            for vidx from uidx+1 to nops(termvars) do
                v := op(vidx, termvars);
                l2 := op(1, v);

                res := res *
                a[lcm(l1, l2)]
                ^((l1*l2/lcm(l1, l2))*
                  degree(term, u)*degree(term, v));
            od;
        od;

        # edges on different cycles of the same size
        for u in termvars do
            l1 := op(1, u); inst1 := degree(term, u);
            # a[l1]^(1/2*inst1*(inst1-1)*l1*l1/l1)
            res := res *
            a[l1]^(1/2*inst1*(inst1-1)*l1);
        od;

        # edges on identical cycles of some size
        for u in termvars do
            l1 := op(1, u); inst1 := degree(term, u);
            if type(l1, odd) then
                # a[l1]^(1/2*l1*(l1-1)/l1);
                res := res *
                (a[l1]^(1/2*(l1-1)))^inst1;
            else
                # a[l1/2]^(l1/2/(l1/2))*a[l1]^(1/2*l1*(l1-2)/l1)
                res := res *
                (a[l1/2]*a[l1]^(1/2*(l1-2)))^inst1;
            fi;
        od;

        all := all + lcoeff(term)*res;
    od;

    all;
end;

kn_pge :=
proc(n,m)
option remember;
local idx_slots, idx_colors, res, term_a, term_b,
    v_a, v_b, inst_a, inst_b, len_a, len_b, p, q;

    if n = 1 or m = 1 then return 1 fi;

    idx_slots := pet_cycleind_edg(n);
    idx_colors := pet_cycleind_symm(m);

    res := 0;

    for term_a in idx_slots do
        for term_b in idx_colors do
            p := 1;

            for v_a in indets(term_a) do
                len_a := op(1, v_a);
                inst_a := degree(term_a, v_a);

                q := 0;

                for v_b in indets(term_b) do
                    len_b := op(1, v_b);
                    inst_b := degree(term_b, v_b);

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

                p := p*q^inst_a;
            od;

            res := res +
            lcoeff(term_a)*lcoeff(term_b)*p;
        od;
    od;

    res;
end;