2019南京网络赛 E 杜教筛+莫比乌斯反演+线性筛

链接： https://nanti.jisuanke.com/t/41302

由于F(p)=p^2-1 且F是积性函数 所以杜教筛的部分其实也可以用min25求

//南京网络赛2019 Ksum https://nanti.jisuanke.com/t/41302
//莫比乌斯反演+线性筛+杜教筛+大数取模+扩展欧拉定理+等比数列求和+数论分块
#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=5e6+5;
const ll mod=1e9+7;
ll inv6=166666668;
ll pri[maxn];bool isp[maxn];int mu[maxn];ll sum[maxn];
ll fuck[maxn];
int p=0;
void init(){//预处理素数跟mu 

    mu[1]=1;isp[1]=1;sum[0]=0;fuck[1]=1;
    for(int i=2;i<maxn;i++){
        if(isp[i]==0)pri[++p]=i,mu[i]=-1,fuck[i]=(1ll*i*i%mod-1+mod)%mod;
        for(int j=1;j<=p&&i*pri[j]<maxn;j++){
            isp[i*pri[j]]=1;
            if(i%pri[j]==0){
                fuck[i*pri[j]]=fuck[i]*pri[j]%mod*pri[j]%mod;
                mu[i*pri[j]]=0;
                break;
            }
            else mu[i*pri[j]]=-mu[i],fuck[i*pri[j]]=1ll*fuck[i]*fuck[pri[j]]%mod;   
        }   
    }   
    for(int i=1;i<maxn;i++)sum[i]=(sum[i-1]+fuck[i])%mod;
} 
map<ll,ll>mp;
ll dfs(ll x){//杜教筛求F的前缀和 
    if(x<=1000000)return sum[x];
    if(mp.count(x))return mp[x];
    ll ans=x*(x+1)%mod*(2*x+1)%mod*inv6%mod;//h的前缀和
    ll nxt=0;
    for(ll i=2;i<=x;i=nxt+1){
        nxt=x/(x/i);
        ans-=1ll*(nxt-i+1)*dfs(x/i)%mod;//要改
        ans=(ans+mod)%mod;
        if(nxt==x)break;
    }
    return mp[x]=ans;
}

ll baoli(ll n){
    ll ans=0;
    for(int i=1;i<=n;i++){
        for(int d=1;d<=i;d++){
            if(i%d==0)ans+=1ll*d*d%mod*mu[i/d]%mod;
        }
    }
    return ans;
}

ll qmod(ll a,ll b){
    ll res=1;
    while(b){
        if(b&1)res=res*a%mod;
        a=a*a%mod;b=b>>1;
    }return res;
}
ll bi(ll x,ll k1,ll k2){//首项为x 公比为x 一共k项

     if(x==1)return (k1-1+mod)%mod;
     
     ll ans=x*(qmod(x,k2)-1+mod)%mod;
     ans=ans*qmod((x-1+mod)%mod,mod-2 )%mod;
     return (ans-x+mod)%mod;
}


int main(){
    //cout<<qmod(6,mod-2)<<endl;
    init();
    /*
    srand(time(0));
    int n=1000+rand()%500;
    cout<<baoli(n)<<" "<<sum[n]<<" "<<dfs(n)<<endl;
    */

    int T;cin>>T;
    while(T--){
        ll n;ll k1=0,k2=0;cin>>n;
        string r;
        getchar();
        cin>>r;getchar();
        int t=1;
        for(int i=0;i<r.size();i++){
            k2=(k2*10+r[i]-'0')%(mod-1);//次方用，取欧拉模
            k1=(k1*10+r[i]-'0')%mod;//平时用，取正常模
        }
        ll ans=0;ll nxt=0;
        for(ll i=1;i<=n;i=nxt+1){
            nxt=n/(n/i);
            ans+=(bi(n/i,k1,k2)*(dfs(nxt)-dfs(i-1)+mod)%mod)%mod;ans%=mod;
        }
        cout<<(ans+mod)%mod<<endl;
    }

    //system("pause");
    return 0;
}