HDU-5780 gcd

题目大意：

求∑gcd(x^a-1, x^b-1)对1e9+7取模的值

解题思路：

官方题解：首先有等式({x}^{a}-1x​a​​−1,{x}^{b}-1x​b​​−1)={x}^{gcd(a,b)}-1x​gcd(a,b)​​−1成立，相当于计算sum sum {x}^{gcd(a,b)}-1∑∑x​gcd(a,b)​​−1 。记s[d]=最大公约数为d的对数，答案sum s[d]*({x}^{d}-1)∑s[d]∗(x​d​​−1). 先用筛法计算出欧拉函数。把正方形分成上三角和下三角计算，最大公约数为d的数量s[d]=2*(phi[1]+phi[2]+...+phi[n/d])-1,对欧拉函数求一个前缀和，直接枚举最大公约数d复杂度为O（n）。观察s[d]计算公式，可以发现对不同的d，若n/d相同，s[d]不发生变化。根据s[d]分段计算，相同的一段的{x}^{d}-1x​d​​−1可以用等比公式求。复杂度为（n+Tsqrt{n}lognn+T√​n​​​logn)

代码：

有个地方需要注意，在求除法取模的时候，需要将除数替换成他对模运算的逆元，否则会出错

#include <cstdio>
#include <cstring>
using namespace std;

typedef long long LL;

const int maxn = 1e6 + 10;
const int mod = 1000000007;

LL cont[maxn];
bool check[maxn];
int tot, prime[maxn], phi[maxn];

void getEuler(){
    memset(check, false, sizeof(check));
    phi[1] = 1;
    tot = 0;
    for(int i = 2; i < maxn; i++){
        if(!check[i]){
            prime[tot++] = i;
            phi[i] = i - 1;
        }
        for(int j = 0; j < tot; j++) {
            if(i * prime[j] > maxn) break;
            check[i * prime[j]] = true;
            if( i % prime[j] == 0){
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }else{
                phi[i * prime[j]] = phi[i] * (prime[j] - 1);
            }
        }
    }

    cont[0] = 0;
    for(int i = 1; i < maxn; ++i){
        cont[i] = (cont[i-1] + phi[i]) % mod;
    }
}
LL Mi(LL a, LL b){
    LL ans = 1;
    while(b){
        if(b & 1) ans = (ans * a) % mod;
        a = (a * a) % mod;
        b >>= 1;
    }
    return ans % mod;
}
LL Cal(LL x, LL be, LL en){
    if(x == 1) return (en - be + 1);
    LL ans = Mi(x, be);
    if(ans < 0) ans += mod;
    LL n = en - be + 1;
    LL tmp = (Mi(x, n) - 1) % mod;
    if(tmp < 0) tmp += mod;
    tmp = (tmp * Mi(x - 1, mod - 2)) % mod;
    ans = (ans * tmp) % mod;
    return ans;
}
int main(){
    int t;
    LL x, n;
    getEuler();
    scanf("%d", &t);
    while(t--){
        scanf("%I64d%I64d", &x, &n);
        if(x == 1) { puts("0");continue; }
        LL nxt, tmp, ans = 0;
        for(int i = 1; i <= n; i = nxt + 1){
            nxt = n / (n / i);
            tmp = (2LL * cont[n / i] - 1) % mod;
            ans = (ans + tmp * Cal(x, i, nxt)) % mod;
            if(ans < 0) ans += mod;
        }
        ans = (ans - n * n) % mod;
        if(ans < 0) ans += mod;
        printf("%I64d
", ans);
    }
    return 0;
}