【题解】[51Nod 1847] 奇怪的数学题【min_25筛 杜教筛 莫比乌斯反演】

题目链接

题意

设 (f(n)) 为 (i) 的第二大因数（(f(n)=dfrac{n}{min(n)})），求 (sum_{i=1}^n sum_{j=1}^n f(gcd(i,j))^k)。(nleq 10^9)，(kleq 50)。模 (2^{32})。

题解

简单莫反可得：

[ans=sum_{i=2}^{n}f(i)^k(-1+2sum_{j=1}^{lfloor frac{n}{i} floor}varphi(i)) ]

整除分块显然，(varphi) 的前缀和可以杜教筛求出，问题在于 (f(i)^k) 的前缀和。

注意到 (f(i)^k) 和 (i^k) 相差的是 (min(i)^k)，而 min_25 筛便将 (min(i)) 不同的数分开来。观察 min_25 筛第一部分的式子（节选）：

[g(n,i-1)-f'(p_i) color{red}left(g(lfloor {nover p_i} floor,i-1)-sum_{j=1}^{i-1}f'(p_j) ight) ]

红色部分便是 (min(j)=p_i) 的合数 (j) 的 (f(j)^k)。

质数要单独拿出来统计。

注意到模数十分阴间，预处理自然数幂和时最好用斯特林数+下降幂。

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int M=7e4+10,L=2e6+10,K=57;
const ll N=1e9+10;
#define uint unsigned int
uint n,k,sqr,lim;
uint pri[L+10],sp1[M],sp2[M],boo[L+10],ppow[M],cnt=0;
uint id1[M],id2[M],g1[M],g2[M],g[M],w[M];
uint fs[M];
uint phi[L+10],phi2[L+10];

int _(ll x){ return x<sqr?id1[x]:id2[n/x]; }
uint qpow(uint x,uint y){
    uint ans=1;
    while(y){
        if(y&1)ans=ans*x;
        x=x*x;
        y>>=1;
    }
    return ans;
}

uint s2[K][K];
void init(int k){
    s2[0][0]=1;
    for(int i=1;i<=k;i++)
        for(int j=1;j<=i;j++)
            s2[i][j]=s2[i-1][j]*j+s2[i-1][j-1];
}
uint s(uint x){
    // sum_{i=0}^{x-1}
    uint ans=0;
    for(int i=0;i<=k;i++){
        uint p=1;
        for(int j=0;j<=i;j++){
            if((x-j)%(i+1)==0)p*=(x-j)/(i+1);
            else p*=x-j;
        }
        ans+=p*s2[k][i];
    }
    return ans;
}
uint sphi(uint x){
    if(x<=lim)return phi[x];
    if(phi2[n/x])return phi2[n/x];
    uint ans=x*1ll*(x+1)/2;
    for(int l=2,r;l<=x;l=r+1){
        r=x/(x/l);
        ans-=(r-l+1)*sphi(x/l);
    }
    // cerr<<"! "<<x<<" "<<ans<<endl;
    return phi2[n/x]=ans;
}
int main(){
    cin>>n>>k;
    init(k);
    sqr=sqrt(n+1);
    for(int i=2;i<=sqr;i++){
        if(!boo[i]){
            pri[++cnt]=i;
            ppow[cnt]=qpow(i,k);
            sp1[cnt]=(sp1[cnt-1]+ppow[cnt]);
            sp2[cnt]=(sp2[cnt-1]+1);
            phi[i]=i-1;
        }
        for(int j=1;j<=cnt&&i*pri[j]<=sqr;j++){
            boo[i*pri[j]]=1;
            if(i%pri[j]==0)break;
        }
    }
    int cnt_=cnt;
    phi[1]=1;
    lim=pow(n+1,0.68);
    cnt=0;
    for(int i=2;i<=lim;i++){
        if(!boo[i]){
            pri[++cnt]=i;
            phi[i]=i-1;
        }
        for(int j=1;j<=cnt&&i*pri[j]<=lim;j++){
            boo[i*pri[j]]=1;
            if(i%pri[j])phi[i*pri[j]]=phi[i]*(pri[j]-1);
            else{
                phi[i*pri[j]]=phi[i]*pri[j];
                break;
            }
        }
    }
    for(int i=1;i<=lim;i++)phi[i]+=phi[i-1];

    int m=1;
    for(ll l=1,r;l<=n;l=r+1,m++){
        r=n/(n/l);
        uint t=n/l;
        w[m]=t;
        g1[m]=s(t+1)-1;
        g2[m]=t-1;
        if(t<sqr)id1[t]=m;
        else id2[n/t]=m;
    }
    --m;
    for(int i=1;i<=cnt_;i++){
        int k=1;
        for(int j=1;j<=m&&w[j]>=pri[i]*1ll*pri[i];j++){
            ll r=_(w[j]/pri[i]);
            fs[j]+=(g1[r]-sp1[i-1]);
            g1[j]-=ppow[i]*(g1[r]-sp1[i-1]);
            g2[j]-=(g2[r]-sp2[i-1]);
        }
    }
    for(int i=1;i<=m;i++)fs[i]+=g2[i];
    uint ans=0;
    for(int i=1;i<=m;i++)
        ans+=(fs[i]-fs[i+1])*(sphi(n/w[i])*2-1);
    cout<<ans;
}