hdu5780 gcd

题意：
已知x和n，求$∑gcd(x^a−1,x^b−1) (1≤a,b≤n)$

思路

记$T(t)=1+x+x^2+...+x^{t-1}$
$x^a-1 = (x-1)(1+x+x^2+...+x^{a-1})=(x-1)*T(a)$
则 $gcd(x^a−1,x^b−1) = (x-1) * gcd(T(a), T(b))$
有公式：$gcd(T(t_1),T(t_2)) = T(gcd(t_1,t_2))$
所以 $gcd(x^a−1,x^b−1) = (x-1) * T(gcd(a,b))=x^{gcd(a,b)}-1$
$∑gcd(x^a−1,x^b−1) =∑(x^{gcd(a,b)}-1)=-n^2+∑x^{gcd(a,b)}$
(a,b)的种数在$n^2$级，但gcd(a,b)的种数在$n$级，假设$gcd(a,b)=i$的对数有$f(i)$个则
$∑x^{gcd(a,b)}（1≤a,b≤n）=∑x^if(n/i)（1≤i≤n）$
$f(i)$很容易由$phi$函数得到，所以整个可以在$O(n)$时间内计算出来，这样T组数据复杂度为$O(Tn)$的，需要优化。
由于$∑x^if(n/i)$中$n/i$的种数在$O(\sqrt{n})$级，所以可以分段计算$∑x^if(n/i)（1≤i≤n）=∑(x^{l_t}+x^{l_t+1}+...+x^{r_t})f(t)（t是n/i取得的所有可能的值）$
对每个t可以$o(1)$算出来，但等比数列求和需要求逆元花费$O(log{n})$时间，所以这一步复杂度为$O(\sqrt(n)log{n})$，加上前面的各种预处理总复杂度为$O(n+T\sqrt{n} log{n})$
对分段方法的说明：
假定当前计算到了i，则需要求出最大的j使得: $\lfloor n/i \rfloor =\lfloor n/j \rfloor$，有一个公式是：$j=\lfloor n/\lfloor n/i \rfloor \rfloor)$ 证明如下
令$n=pi+q$，则
$\lfloor n/i \rfloor=\lfloor n/j\rfloor =p$ ———A
$\lfloor n/(j+1) \rfloor <p$ ———B
由A得$\lfloor(pi+q)/j\rfloor =p$，所以$pi+q≥pj$ ———-C
由B得$\lfloor(pi+q)/(j+1) \rfloor <p$，所以$pi+q<p(j+1)$ ———D
综合C、D得 $i+q/p-1<j≤i+q/p$ 由于左右边界刚好相差1，所以区间内刚好有1个整数，即$j=i+\lfloor q/p\rfloor$
$i+\lfloor q/p\rfloor = i+\lfloor (n-pi)/p\rfloor = i+\lfloor n/p\rfloor -i=\lfloor n/p\rfloor=\lfloor n/\lfloor n/i\rfloor\rfloor$

#include <bits/stdc++.h>
using namespace std;
const int mod = 1e9 + 7;
const int N = 1e6 + 7;
typedef long long ll;
int phi[N], f[N], g[N], x, n, _, sqrtn, ans, l, r;

void init() {
    for (int i = 1; i < N; i++) phi[i] = i;
    for (int i = 2; i < N; i++) {
        if (phi[i] == i) {
            for (int j = i; j < N; j += i) {
                phi[j] = phi[j] / i * (i - 1);
            }
        }
    }
    for (int i = 1; i < N; i++) {
        f[i] = (f[i - 1] + phi[i]) % mod;
        g[i] = (2ll * f[i] - 1 + mod) % mod;
    }
}
ll powermod(ll a, ll n, ll mod) {
    ll ans = 1;
    while (n) {
        if (n & 1) ans = ans * a % mod;
        a = a * a % mod;
        n >>= 1;
    }
    return ans;
}
ll F(ll l, ll r) {
    if (x == 1) return r - l + 1;
    return (powermod(x, r + 1, mod) - powermod(x, l, mod) + mod) % mod * powermod(x - 1, mod - 2, mod) % mod;
}
int main() {
    init();
    cin >> _;
    while (_--) {
        cin >> x >> n;
        ans = 0;
        l = r = 1;
        for (; l <= n; l = r + 1) {
            r = n / (n / l);
            ans = (ans + g[n / l] * F(l, r) % mod) % mod;
        }
        cout << (ans - (ll)n * n % mod + mod) % mod << endl;
    }
    return 0;
}