题目大意:
求∑gcd(x^a-1, x^b-1)对1e9+7取模的值
解题思路:
官方题解:首先有等式(xa−1,xb−1)=xgcd(a,b)−1成立,相当于计算∑∑xgcd(a,b)−1 。记s[d]=最大公约数为d的对数,答案∑s[d]∗(xd−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]分段计算,相同的一段的xd−1可以用等比公式求。复杂度为(n+T√nlogn)
代码:
有个地方需要注意,在求除法取模的时候,需要将除数替换成他对模运算的逆元,否则会出错
#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;
}