思路
很明显的一个思路：先搞出有多少种珠子，再求有多少种项链。
珠子
考虑这个式子：
\[
S3=\sum_{i=1}^a \sum_{j=1}^a\sum_{k=1}^a [\gcd(i,j,k)==1]
\]
显然可以莫比乌斯反演一波，但这个是对的吗？
当有两个数字相同时只被算了3遍，而三个都相同的只被算了一遍。
\[
S2=\sum_{i=1}^a\sum_{j=1}^a [\gcd(i,j)==1]
\]
显然有\(S1=1\)，那么就会得到最终答案：
\[
ans=\frac 1 6 (S3+3S2+2)
\]
项链
既然要求旋转之后不相同，那么自然想到polya定理。
设不同珠子的个数为\(m\)，那么可以得到
\[
ans=\frac 1 n \sum_{d|n}f(d)\varphi(\frac n d)
\]
其中\(f(d)\)表示\(d\)个珠子串成项链，使得相邻的不相等的概率，注意不能旋转，且\(f(1)=0\)。
考虑从右往左数第二个是否和从左往右数第一个相同，那么有
\[
f_n=(m-2)f_{n-1}+(m-1)f_{n-2}
\]
边界条件：\(f_1=0,f_2=m(m-1)\)。
可以看出这是个线性齐次XXXX递推式，用特征根可以得到
\[
x_1=-1,x_2=m-1\\
\alpha = m-1,\beta=1
\]
所以得到\(f_d=(m-1)^d+(-1)^d(m-1)\)，可以\(O(\log d)\)得到了。
于是\(ans\)似乎也可以很快得到了。
汇总
总结一下算法流程：先整除分块处理出\(m\)，然后polya定理搞出答案。
有一个问题：\(\varphi(\frac n d)\)怎么快速算出来？
看大佬的做法，都是先给\(n\)分解质因数，然后\(dfs\)枚举每种组合，一次性把\(\varphi(d)\)全都求出来，复杂度似乎是\(O(n的因数个数)\)，也就不超过\(O(\sqrt n)\)。
还有一个问题：\(n\)是模数的倍数时怎么办？
改成\(MOD=mod^2\)，就不会出现这种情况了，最后输出答案时好像还要用奇奇怪怪的方法搞一搞。
代码
咕咕咕咕咕咕
以后再写吧。
