ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

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

2020-12-24 23:01:22  阅读:159  来源: 互联网

标签:10 min int 题解 sum uint 25 ll


题目链接

题意

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

题解

简单莫反可得:

\[ans=\sum_{i=2}^{n}f(i)^k(-1+2\sum_{j=1}^{\lfloor \frac{n}{i}\rfloor}\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 {n\over p_i}\rfloor,i-1)-\sum_{j=1}^{i-1}f'(p_j) \right) \]

红色部分便是 \(\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;
}

标签:10,min,int,题解,sum,uint,25,ll
来源: https://www.cnblogs.com/wallbreaker5th/p/14186813.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有