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+T\sqrt{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\n", ans);
}
return 0;
}