1. def countNecklaces(n, m):
2.     N = n*m
3.     R = PolynomialRing(QQ, m, 't')
4.     ts = R.gens()
5.     f = sum(t for t in ts)
6.     R2 = PolynomialRing(R, N, 'a')
7.     aGens = R2.gens()
8.     ds = divisors(N)
9.     F = 1/N * sum(euler_phi(N/d) * aGens[N/d-1]**d for d in ds)
10.     g = F(*[f(*[t**k for t in ts]) for k in range(1, N+1)])
11.     #print (g)
12.     wantMono = prod(t**n for t in ts)
13.     return g.coefficient(wantMono)
14.  
15.  
16. def countBracelets(n, m):
17.     N = n*m
18.     R = PolynomialRing(QQ, m, 't')
19.     ts = R.gens()
20.     f = sum(t for t in ts)
21.     R2 = PolynomialRing(R, N, 'a')
22.     aGens = R2.gens()
23.     ds = divisors(N)
24.     F = (1/(2*N) * sum(euler_phi(N/d) * aGens[N/d-1]**d for d in ds)
25.          + ( (1/2*aGens[0]*aGens[1]**int((N-1)/2)) if N%2==1
26.              else (1/4*(aGens[0]**2 * aGens[1]**int((N-2)/2) + aGens[1]**int(N/2)))
27.            )
28.         )
29.    
30.     g = F(*[f(*[t**k for t in ts]) for k in range(1, N+1)])
31.     #print (g)
32.     wantMono = prod(t**n for t in ts)
33.     return g.coefficient(wantMono)
34.  
35.  
36. print (countNecklaces(20, 3))
37. print (countBracelets(20, 3))
38.  
39. #print ([countBracelets(k, 3) for k in range(1, 10)])