with(numtheory);

Y :=
proc(n, k)
option remember;
local d, dd, ind, orbit, orbits, pos, shft;

    if k = 1 then return 1 fi;

    orbits := {};

    for ind from k^n to 2*k^n-1 do
        d := convert(ind, base, k);

        dd := [seq(d[q], q=1..n), d[1], d[2]];

        for pos to n do
            if dd[pos] = 1 and
            dd[pos+1] = 1 and dd[pos+2] = 0 then
                break;
            fi;
        od;

        if pos = n+1 then
            orbit := {};

            for shft to n do
                orbit :=
                {op(orbit),
                 [seq(d[p], p=shft .. n),
                  seq(d[p], p=1..shft-1)]};
            od;

            orbits := {op(orbits), orbit};
        fi;
    od;

    nops(orbits);
end;


pet_varinto_cind :=
proc(poly, ind)
local subs1, subs2, polyvars, indvars, v, pot, res;

    res := ind;

    polyvars := indets(poly);
    indvars := indets(ind);

    for v in indvars do
        pot := op(1, v);

        subs1 :=
        [seq(polyvars[k]=polyvars[k]^pot,
             k=1..nops(polyvars))];

        subs2 := [v=subs(subs1, poly)];

        res := subs(subs2, res);
    od;

    res;
end;

pet_cycleind_cyclic :=
proc(n)
local d, s;

    s := 0;
    for d in divisors(n) do
        s := s + phi(d)*a[d]^(n/d);
    od;

    s/n;
end;

pet_cycleind_cyclic_num :=
proc(n, k)
local d, s;

    s := 0;
    for d in divisors(n) do
        s := s + phi(d)*k^(n/d);
    od;

    s/n;
end;


R :=
proc(n, k)
option remember;
local res, seg, inv, ovars, ofree, rest, gf;

    if k = 1 then return 1 fi;

    ovars := add(Z*Y[q], q=0..k-2);
    ofree := add(Y[q], q=0..k-2);

    inv := add(W*Z*ovars^q, q=1..n)
    + add(add(W^p*Z^p*(ovars-Z*Y[0])*ovars^q, q=0..n),
          p=2..n);

    res := 0;

    for seg to floor(n/2) do
        gf := expand(pet_varinto_cind
                     (inv, pet_cycleind_cyclic(seg)));
        res := res + coeff(gf, Z, n);
    od;

    rest := pet_varinto_cind(ofree, pet_cycleind_cyclic(n));

    1 + subs({W=1, seq(Y[q]=1, q=0..k-2)}, res + rest);
end;

R2 :=
proc(n, k)
option remember;
local seg, gf, res;

    if k = 1 then return 1 fi;

    res := 0;

    gf := z*add((k-1)^q*z^q, q=1..n)
    + add(z^p*(k-2)*z*add((k-1)^q*z^q, q=0..n), p=2..n);

    for seg to floor(n/2) do
        res := res +
        coeff(expand(pet_varinto_cind
                     (gf, pet_cycleind_cyclic(seg))),
              z, n);
    od;


    1 + pet_cycleind_cyclic_num(n, k-1) + res;
end;