def countNecklaces(n, m):
    N = n*m
    R = PolynomialRing(QQ, m, 't')
    ts = R.gens()
    f = sum(t for t in ts)
    R2 = PolynomialRing(R, N, 'a')
    aGens = R2.gens()
    ds = divisors(N)
    F = 1/N * sum(euler_phi(N/d) * aGens[N/d-1]**d for d in ds)
    g = F(*[f(*[t**k for t in ts]) for k in range(1, N+1)])
    #print (g)
    wantMono = prod(t**n for t in ts)
    return g.coefficient(wantMono)


sol = countNecklaces(20, 3)
print (sol)