題目傳送門：洛谷P3307。這題在bzoj上是許可權題。

題意簡述：

這題分為兩個部分：

① 有一些珠子，每個珠子可以看成一個無序三元組。三元組要滿足三個數都在$1$到$m$之間，並且三個數互質，兩個珠子不同當且僅當這個三元組不同。計算有多少種不同的珠子。

② 把這些珠子串成一個環，要滿足相鄰的珠子不同。兩個環不同當且僅當旋轉任意角度後仍然不同。計算有多少種不同的環。

題解：

分成兩部分做。

第一部分：

考慮計算三元組的個數，轉無序為有序，再去重。

答案=(三個都不同的有序三元組方案)/6+(兩個相同，另一個不同的方案)/3+(三個都相同的方案)。

容斥一下得到答案=(三元組的方案+二元組的方案*3+一元組的方案*2)/6。

因為一元組只有(1)滿足條件，所以答案是(2+三元組的方案+二元組的方案*3)/6。

考慮如何求出兩種方案。

三元組的方案是\(\sum_{i=1}^m\sum_{j=1}^m\sum_{k=1}^m[\gcd(i,j,k)=1]\)，二元組同理。

顯然是莫反套路，三元組的答案是\(\sum_{d=1}^m\mu(d){\left\lfloor\frac{m}{d}\right\rfloor}^3\)，二元組同理。

數論分塊求出答案即可，最後乘上6的逆元。這一步複雜度$\Theta(m+T\sqrt{m})$。

第二部分：

知道了不同珠子的數量，要求出本質不同的環的個數。

Burnside引理套路。最終方案數等於每個置換的不動點個數的平均數，即\(\frac{1}{n}\sum_{i=1}^nf(i)\)，\(f(i)\)表示旋轉\(i\)格的不動點數量。

稍微化簡一下：\(\frac{1}{n}\sum_{d|n}\varphi(\frac{n}{d})f(d)\)。

考慮計算\(f(x)\)，當$x$是$n$的因數時，$f(x)$就等於不考慮旋轉時的長度為$x$的環的數量。

假設不同珠子的數量為\(k\)，不加證明地給出一個式子：\(f(x)=(k-1)^x+(-1)^x(k-1)\)。這個式子可以遞推得出。

那麼根據這個式子和上面的式子計算即可。

要注意\(n\)太大了，要求出\(\varphi\)的值比較困難，考慮DFS它的每個質因數，按照\(\varphi\)是個積性函式以及公式，求得\(\varphi\)。

要注意，最後除掉\(n\)的時候，\(n\)可能是模數的倍數導致沒有逆元。可以發現\(n\)不會是模數平方的倍數，所以把模數平方後再做一遍，最後除掉模數這個因子即可。

  1 #include <cstdio>
  2 
  3 #define reg register
  4 typedef unsigned long long ULL;
  5 const ULL MOD = 1000000007ll;
  6 const ULL Inv61 = 166666668ll;
  7 const ULL Inv62 = 833333345000000041ll;
  8 ULL Mod;
  9 ULL Inv6;
 10 const int MN = 10000001;
 11 
 12 ULL TN[11];
 13 int TA[11], MA;
 14 
 15 bool ip[MN];
 16 int p[MN], pc;
 17 int mu[MN];
 18 inline void SieveInit() {
 19     ip[0] = ip[1] = 1;
 20     mu[1] = 1;
 21     for (reg int i = 2; i <= MA; ++i) {
 22         if (!ip[i])
 23             p[++pc] = i,
 24             mu[i] = -1;
 25         for (reg int j = 1; j <= pc; ++j) {
 26             reg int k = p[j] * i;
 27             if (k > MA) break;
 28             ip[k] = 1;
 29             if (i % p[j]) mu[k] = -mu[i];
 30             else break;
 31         }
 32     }
 33     for (reg int i = 2; i <= MA; ++i)
 34         mu[i] += mu[i - 1];
 35 }
 36 
 37 int O;
 38 inline ULL Mul(ULL x, ULL y) {
 39     if (!O) return x * y % Mod;
 40     return (x * y - (ULL)((long double) x / Mod * y) * Mod + Mod) % Mod;
 41 }
 42 
 43 ULL N; int A;
 44 ULL M;
 45 inline void SolveM() {
 46     M = 2;
 47     for (reg int i = 1, j, k; i <= A; i = j + 1) {
 48         k = A / i, j = A / k;
 49         M = (M + Mul(Mul(Mul(k, k), k + 3), (mu[j] - mu[i - 1] + Mod) % Mod)) % Mod;
 50     }
 51     M = Mul(M, Inv6);
 52 }
 53 
 54 ULL Pow[60];
 55 inline void PowInit() {
 56     Pow[0] = M - 1;
 57     for (reg int i = 1; i < 60; ++i) Pow[i] = Mul(Pow[i - 1], Pow[i - 1]);
 58 }
 59 inline ULL qPow(ULL E) {
 60     ULL A = 1;
 61     for (reg int j = 0; E; E >>= 1, ++j)
 62         if (E & 1) A = Mul(A, Pow[j]);
 63     return A;
 64 }
 65 inline ULL Inv(ULL B) {
 66     ULL A = 1;
 67     for (reg ULL E = MOD - 2; E; E >>= 1, B = B * B % MOD)
 68         if (E & 1) A = A * B % MOD;
 69     return A;
 70 }
 71 
 72 ULL b[15]; int e[15], cnt;
 73 ULL Ans;
 74 inline ULL F(ULL x) {
 75     return (qPow(x) + (x & 1 ? Mod - M + 1 : M - 1)) % Mod;
 76 }
 77 void DFS(int st, ULL now, ULL phi) {
 78     if (st > cnt) {
 79         Ans = (Ans + Mul(phi % Mod, F(N / now))) % Mod;
 80         return;
 81     }
 82     DFS(st + 1, now, phi);
 83     for (reg int i = 1; i <= e[st]; ++i) {
 84         now *= b[st];
 85         phi *= i == 1 ? b[st] - 1 : b[st];
 86         DFS(st + 1, now, phi);
 87     }
 88 }
 89 inline ULL Solve() {
 90     ULL NN = N; cnt = 0;
 91     for (reg ULL i = 2; i * i <= NN; ++i) if (NN % i == 0) {
 92         b[++cnt] = i, e[cnt] = 0;
 93         while (NN % i == 0) NN /= i, ++e[cnt];
 94     } if (NN > 1) b[++cnt] = NN, e[cnt] = 1;
 95     Ans = 0; DFS(1, 1, 1);
 96     if (O) Ans = Ans / MOD * Inv(N / MOD) % MOD;
 97     else Ans = Ans * Inv(N % MOD) % MOD;
 98     return Ans;
 99 }
100 
101 int main() {
102     int Tests;
103     scanf("%d", &Tests);
104     for (int i = 1; i <= Tests; ++i)
105         scanf("%llu%d", TN + i, TA + i),
106         MA = TA[i] > MA ? TA[i] : MA;
107     SieveInit();
108     for (int i = 1; i <= Tests; ++i) {
109         N = TN[i], A = TA[i];
110         O = N % MOD ? 0 : 1;
111         if (O) Mod = MOD * MOD, Inv6 = Inv62;
112         else Mod = MOD, Inv6 = Inv61;
113         SolveM();
114         PowInit();
115         printf("%llu\n", Solve());
116     }
117     return 0;
118 }
119 
120 // 1. 求出本質不同的珠子數量，容斥 + 莫比烏斯反演 + 數論分塊
121 // 2. 求出答案，Burnside 引理 + 數論分塊