[BZOJ4407]于神之怒加强版 CODE

#include<bits/stdc++.h>
#define ll long long
#define For(i,a,b) for(ll i=(a);i<=(b);++i)
#define Rep(i,a,b) for(ll i=(a);i>=(b);--i)
const ll N=1e6+10;
using namespace std;
const ll mod=1e9+7;

ll vis[N],tot,p[N];
void init(ll n){//质数筛
    For(i,2,n){
        if(!vis[i])p[++tot]=i;
        For(j,1,tot){
            if(i*p[j]>n)break;
            vis[i*p[j]]=1;
            if(i%p[j]==0)break;
        }
    }
}

ll cnt,t[N],num[N];//t为质因数，num为出现的个数
void solve(ll x){//质因数分解
    for(ll i=1;p[i]*p[i]<=x;i++){
        if(x%p[i]==0){
            t[++cnt]=p[i];
            num[cnt]=0;
            while(x%p[i]==0)num[cnt]++,x/=p[i];
        }
    }
    if(x>1)t[++cnt]=x,num[cnt]=1;
}

ll a,b,ans;
void dfs(ll x,ll d,ll s){//第x个质因数，目前枚举的因数为d，f(d)=s
    if(x>cnt){//按公式计算答案
        ll t1=b/d,t2=(a+d-1)/d;
        ans=(ans+(t1+t2)*(t1-t2+1)*s%mod+mod)%mod;
        return;
    }
    dfs(x+1,d,s);//不选当前质因数
    s=s*(1-t[x]);//积性函数直接计算
    For(i,1,num[x]){//选i个当前质因数
        d*=t[x];
        dfs(x+1,d,s);
    }
}

void mian(){

    scanf("%lld%lld",&a,&b);
    cnt=ans=0;//记得清空
    solve(b);
    dfs(1,1,1);
    printf("%lld\n",ans*500000004%mod*b%mod);//500000004是2在模1000000007下的逆元

}

int main(){
    init(50000);
    int T=1;
    scanf("%d",&T);
    while(T--)mian();
    return 0;
}