pet_cycleind_symm :=
proc(n)
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_cycleind_knn :=
proc(n)
option remember;
local cindA, cindB, sind, t1, t2, srcterm, term,
    res, len, len1, len2, cycs, uidx, vidx, interlv,
    u, v, inst, q, v1, v2;

    if n=1 then
        sind := [a[1]];
    else
        sind := pet_cycleind_symm(n);
    fi;

    cindA := 0;

    for t1 in sind do
        for t2 in sind do
            q := 1;

            for v1 in indets(t1) do
                len1 := op(1, v1);

                for v2 in indets(t2) do
                    len2 := op(1, v2);

                    len := lcm(len1, len2);
                    q := q *
                    a[len]^((len1*len2/len) *
                    degree(t1, v1)*degree(t2, v2));
                od;
            od;

            cindA := cindA +
            lcoeff(t1)*lcoeff(t2)*q;
        od;
    od;

    cindB := 0;

    interlv := {seq(a[q]=a[2*q], q=1..n)};
    for srcterm in sind do
        term := subs(interlv, srcterm);
        termvars := indets(term); res := 1;

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

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

                res := res *
                a[lcm(len1, len2)]
                ^((len1*len2/2/lcm(len1, len2))*
                  degree(term, u)*degree(term, v));
            od;
        od;

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

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

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

    (cindA+cindB)/2;
end;


A007139 :=
proc(n)
    option remember;
    local cind;

    cind := pet_cycleind_knn(n);
    subs([seq(a[p]=2, p=1..n*n)], cind);
end;