题解：HDU5780 gcd

题意：求：

\[\sum_{i=1}^n\sum_{j=1}^n\gcd(x^i-1,x^j-1) \]

\(T\) 次询问，\(T\le300\)，\(x,n\le10^6\)。对 \(10^9+7\) 取模。

---

记录一下这个神奇的东西：

\[\gcd(x^a-1,x^b-1)=x^{\gcd(x,y)}-1 \]

有类似性质的还有斐波那契数列等。

证明：

定义分圆多项式：

\[\phi_n(x)=\prod_{i=0}^{n-1}(x-\omega_n^i)^{[i\perp n]} \]

即以所有 \(n\) 次本原单位根为根的多项式，其次数为 \(\varphi(n)\)，系数为 \(\pm1\)。

有什么用？由于 \(x^n-1\) 的根就是所有 \(n\) 次单位根，有：

\[\begin{aligned} x^n-1&=\prod_{i=0}^{n-1}(x-\omega_n^i)\\ &=\prod_{d|n}\prod_{i=0}^{n-1}(x-\omega_n^i)^{[\gcd(n,i)=d]}\\ &=\prod_{d|n}\prod_{i=0}^{n/d-1}(x-\omega_{n/d}^i)^{[i\perp n/d]}\\ &=\prod_{d|n}\phi_{n/d}(x)\\ &=\prod_{d|n}\phi_d(x) \end{aligned} \]

即 \(x^n-1\) 可以表示成其所有因数的分圆多项式的乘积。

得到了 \(x^n-1\) 的分解后，上面的性质就一目了然了。

---

说回本题，所求即为：

\[\sum_{i=1}^n\sum_{j=1}^n\gcd(x^i-1,x^j-1)=\sum_{i=1}^n\sum_{j=1}^n(x^{\gcd(i,j)}-1) \]

枚举 \(\gcd\)：

\[\sum_{d=1}^n(x^d-1)\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=d] \]

后面那一坨典中典了，就是 \(2\sum\limits_{i=1}^{n/d}\varphi(i)-1\)。

那么原式即为：

\[\sum_{d=1}^n(x^d-1)\left(2\sum_{i=1}^{n/d}\varphi(i)-1\right) \]

整除分块，现在要求：

\[\sum_{d=l}^r(x^d-1)=\dfrac{x^l(x^{r-l+1}-1)}{x-1}-(r-l+1) \]

然后就能做了，时间复杂度为 \(\mathcal O(T\sqrt n\log n)\)，如果预处理 \(x\) 的块幂的话可以做到 \(\mathcal O(T\sqrt n)\)。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e6 + 5, p = 1e9 + 7;
int T, x, n, inv, tot, pri[N], phi[N];
bool vis[N];
void sieve(int n) {phi[1] = 1;for (int i = 2, k, p; i <= n; i++) {if (!vis[i]) pri[++tot] = i, phi[i] = i - 1;for (int j = 1; j <= tot && (k = i * (p = pri[j])) <= n; j++) {vis[k] = 1;if (!(i % p)) { phi[k] = phi[i] * p; break; }phi[k] = phi[i] * (p - 1);}}for (int i = 1; i <= n; i++) phi[i] = (phi[i - 1] + phi[i]) % p;
}
int qpow(int a, int b) {int c = 1;while (b) { if (b & 1) c = (ll)c * a % p; a = (ll)a * a % p, b >>= 1; }return c;
}
int sum(int l, int r) {return ((ll)qpow(x, l) * (qpow(x, r - l + 1) - 1) % p * inv % p + p - (r - l + 1)) % p;
}
int main() {sieve(1e6);scanf("%d", &T);while (T--) {scanf("%d%d", &x, &n);if (x == 1) { printf("0\n"); continue; }inv = qpow(x - 1, p - 2);int ans = 0;for (int l = 1, r, k; l <= n; l = r + 1) {r = n / (k = n / l);ans = (ans + (ll)(phi[k] * 2 - 1) * sum(l, r)) % p;}printf("%d\n", ans);}
}
/*
5
3 1
4 2
8 7
10 5
10 8
*/