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)


def countBracelets(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/(2*N) * sum(euler_phi(N/d) * aGens[N/d-1]**d for d in ds)
         + ( (1/2*aGens[0]*aGens[1]**int((N-1)/2)) if N%2==1
             else (1/4*(aGens[0]**2 * aGens[1]**int((N-2)/2) + aGens[1]**int(N/2)))
           )
        )
    
    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)


print (countNecklaces(20, 3))
print (countBracelets(20, 3))

#print ([countBracelets(k, 3) for k in range(1, 10)])