题意:
给定1<=a<=b<=109,求∑i=ablcm(i,b)?
分析:
∑i=ablcm(i,b)=∑i=abi∗bgcd(i,b)=b∗∑d∑a<=i<=bid∗(gcd(i,b)=d)=b∗∑d|b∑⌈ad⌉<=i<=bdi∗(gcd(i,bd)=1)=b∗∑d|b∑⌈ad⌉<=i<=bdi∗∑d′|gcd(i,bd)μ(d′)=b∗∑d|b∑d′|bdμ(d′)∑⌈ad⌉<=i<=bd,d′|ii=b∗∑d|b∑d′|bdμ(d′)∗d′∗(⌈⌈ad⌉d′⌉+⌊⌊bd⌋d′⌋)∗(⌊⌊bd⌋d′⌋−⌈⌈ad⌉d′⌉+1)2=b∗∑d|b∑d′|bdμ(d′)∗d′∗(⌈add′⌉+⌊bdd′⌋)∗(⌊bdd′⌋−⌈add′⌉+1)2=b∗∑d|b(⌈ad⌉+⌊bd⌋)∗(⌊bd⌋−⌈ad⌉+1)2∑d′|dμ(d′)∗d′
最后一步转化是另 d=dd′ 。
做到最后发现需要枚举
b
的所有约数,此处我用
易知:
f(d)
是一个积性函数,
f(p)=1−p,f(pk)=f(p)
此处
p
为素数。
所以枚举约数
算一下复杂度:
n
分解质因数复杂度不超过
所以总的最大复杂度为
o(T∗(n√+256∗8))
,注意所有数据都达到不了这个复杂度。
下面是代码:
#include <bits/stdc++.h>
#define LL long long
#define FOR(i,x,y) for(int i = x;i < y;++ i)
#define IFOR(i,x,y) for(int i = x;i > y;-- i)
using namespace std;
const LL Mod = 1000000007;
const int maxn = 100010;
const LL inv = 500000004;
LL a,b;
int prime[maxn];
bool check[maxn];
int pri[35],pri_cnt,len[35];
vector <int> fac;
void Mobius(){
memset(check,false,sizeof(check));
prime[0] = 0;
FOR(i,2,maxn){
if(!check[i]){
prime[++prime[0]] = i;
}
FOR(j,1,prime[0]+1){
if(i*prime[j] >= maxn) break;
check[i*prime[j]] = true;
if(i%prime[j] == 0) break;
}
}
}
void dfs(int res,int l){
if(l >= pri_cnt) {fac.push_back(res);return;}
int tem = 1;
dfs(res,l+1);
FOR(i,1,len[l]+1){
tem *= pri[l];
dfs(res*tem,l+1);
}
}
void Get_Fac(){
pri_cnt = 0;
int b_c = b;
FOR(i,1,prime[0]+1){
if(prime[i]*prime[i] > b_c) break;
if(b_c % prime[i] == 0) {pri[pri_cnt++] = prime[i]; len[pri_cnt-1] = 0;}
while(b_c % prime[i] == 0) {b_c /= prime[i];len[pri_cnt-1] ++;}
}
if(b_c > 1) {pri[pri_cnt++] = b_c; len[pri_cnt-1] = 1;}
fac.clear();
dfs(1,0);
}
void work(){
Get_Fac();
LL ans = 0;
FOR(i,0,(int)fac.size()){
int v = fac[i];
LL tt = a+v-1;
LL tem = (((tt/v+b/v)%Mod)*((b/v-tt/v+1+Mod)%Mod)%Mod)*inv%Mod;
LL res = 1;
FOR(j,0,pri_cnt){
if(v % pri[j]) continue;
LL t = (1-pri[j]+Mod)%Mod;
res = (res*t)%Mod;
}
ans += (res*tem)%Mod;
ans %= Mod;
}
ans = (ans*b)%Mod;
printf("%lld\n",ans);
}
int main()
{
//freopen("test.in","r",stdin);
//freopen("out.txt","w",stdout);
Mobius();
int T; scanf("%d",&T);
while(T--){
scanf("%lld%lld",&a,&b);
work();
}
return 0;
}