zoukankan      html  css  js  c++  java
  • 【数论】Min_25筛

    Min_25筛

    学习资料:OI Wiki

    P5325题解 by wucstdio

    用法:

    求解 (sum_{i=1}^n f(i)) ,满足 (f(x)) 是一个积性函数,且 (f(p^e)) 是关于 (p) 的低阶多项式。

    板子&&例题

    例题1: P5325 【模板】Min_25筛

    题意:

    定义积性函数 (f(x)) ,且 (f(p^k)=p^k(p^k-1))(p) 是一个质数)。对给定整数 (n)(1le nle10^{10}) ,求 (sum_{i=1}^nf(i)) 。对 (10^9+7) 取模。

    解:

    首先,对 (i) 按照质数和合数分类:

    [sum^n_{i=1}f(i)=sum_{1le ple n}f(p)+sum^n_{i=1}f(i)[i otin prime] ]

    然后,枚举后边合数的最小质因子以及最小质因子次数。所有合数的最小质因子一定小于等于 (sqrt{n})

    [sum_{1le ple n}f(p)+sum_{1le p^ele n,1le p,le{sqrt{n}}}f(p^e)(sum_{i=1}^{frac n{p^e}}f(i)*[minp>p]) ]

    其中, (minp) 表示 (i) 的最小质因子。

    将式子分为前后两部分。

    (clubsuit) 第一部分:

    设完全积性函数 (g(x))其在质数处取值 (g(x)=f(x))

    设数组 (G(n,j)) 满足:

    [G(n,j)=sum^n_{i=1} g(i)[iin primevee minp>p_j] ]

    因为 (G(n,j)) 可以由 (G(n,j-1)) 减去 最小质因子为 (p_j) 的合数函数和 得到,而最小质因子为 (p_j) 的合数函数和,提取出质因子 (p_j) 后可表示为 :

    [g(j)[G(frac n{p_j},j)-G(p_{j-1},j-1)] ]

    所以有:

    [G(n,j)=G(n,j-1)-g(j)[(G(frac n{p_j},j-1)-G(p_{j-1},j-1))] ]

    其中:

    [G(p_{j-1},j-1)=sum_{i=1}^{j-1}g(p_j)=sum_{i=1}^{j-1}f(p_j) ]

    所以:

    [G(n,j)=G(n,j-1)-g(p_j)[(G(frac n{p_j},j-1)-sum_{i=1}^{j-1}g(p_i)]\ ]

    对于 (G(n,j)) 函数,(1)(n) 中所有质数的函数值和 等值于 (G(n,x)) ,其中 (p_x) 是最后一个小于等于 (sqrt{n}) 的质数。为了方便,计其为 (G(n)) 。可用整除分块求。


    (clubsuit) 第二部分:

    (S(n,x)) 表示求 (1)(n) 中所有最小质因子大于 (p_x) 的函数值之和。则题目的答案就是 (S(n,0))(p_0=1) ) 。

    (S(n,x)) 分成两部分,第一部分是大于 (p_x) 的质数,即 $G(n)-sum_{i=1}^x g(p_i) $,第二部分是最小质因子大于 (p_x) 的合数,枚举最小质因子:

    [S(n,x)=G(n)-sum_{i=1}^x g(p_i)+sum_{p_k^ele nwedge k>x}f(p_k^e)(S(frac n{p_k^e},k)+[e e 1]) ]

    (clubsuit)最后:

    将第一部分定义的积性函数 (g(x)) ,根据 (g(p)=f(p)=p*(p-1)) 的定义,拆成 (g1(x)=x)(g2(x)=x^2) 进行求解。

    #include<bits/stdc++.h>
    #define mem(a,b) memset(a,b,sizeof(a))
    using namespace std;
    typedef long long ll;
    const int maxn=1e6+50;
    const int mod=1e9+7;
    
    const int inv6=166666668;
    const int inv2=500000004;
    int pri[maxn],ncnt,T,m,id1[maxn],id2[maxn];
    bool isp[maxn];
    ll n,a[maxn],g1[maxn],g2[maxn],sum1[maxn],sum2[maxn];
    void init()
    {
        T=sqrt(n+0.5);
        for(int i=2;i<=T;i++){
            if(!isp[i])pri[++ncnt]=i,sum1[ncnt]=(sum1[ncnt-1]+i)%mod,sum2[ncnt]=(sum2[ncnt-1]+1ll*i*i%mod)%mod;
            for(int j=1;j<=ncnt&&i*pri[j]<=T;j++){
                isp[i*pri[j]]=true;
                if(i%pri[j]==0)break;
            }
        }
        for(ll l=1,r;l<=n;l=r+1){
            r=n/(n/l);
            a[++m]=n/l;
            if(a[m]<=T)id1[a[m]]=m;else id2[n/a[m]]=m;
            g1[m]=a[m]%mod;
            g2[m]=(g1[m]*(g1[m]+1)%mod*(g1[m]*2+1)%mod*inv6%mod-1+mod)%mod;
            g1[m]=(g1[m]*(g1[m]+1)%mod*inv2%mod-1+mod)%mod;
        }
        for(int i=1;i<=ncnt;i++){
            ll p2=1ll*pri[i]*pri[i],v;
            for(int j=1,id;j<=m&&p2<=a[j];j++){
                v=a[j]/pri[i];
                id=v<=T?id1[v]:id2[n/v];
                g1[j]=(g1[j]-1ll*pri[i]*(g1[id]-sum1[i-1]+mod)%mod+mod)%mod;
                g2[j]=(g2[j]-1ll*pri[i]*pri[i]%mod*(g2[id]-sum2[i-1]+mod)%mod+mod)%mod;
            }
        }
    }
    ll S(ll x,int y)
    {
        if(pri[y]>=x)return 0;
        int id=x<=T?id1[x]:id2[n/x];
        ll sum=(g2[id]-g1[id]-(sum2[y]-sum1[y]+mod)%mod+mod+mod)%mod;
        for(int k=y+1;k<=ncnt&&1ll*pri[k]*pri[k]<=x;k++){
            ll pe=pri[k];
            for(int e=1;pe<=x;pe*=pri[k],e++)
                sum=(sum+pe%mod*((pe-1)%mod)%mod*(S(x/pe,k)+(e!=1))%mod)%mod;
        }
        return sum;
    }
    int main()
    {
        scanf("%lld",&n);
        init();
        printf("%lld
    ",(S(n,0)+1)%mod);
    }
    



    例题2:#6053.简单的函数

    题意:

    (f(x)) 有如下性质:

    • (f(1)=1)
    • (f(p^c)=pigoplus c)(p) 是质数,(igoplus) 表示异或)
    • (f(ab)=f(a)f(b))(a)(b) 互质)

    对于给定 (n)(1le nle 10^{10}) ,求 (sum^n_{i=1}f(i)) 。对 (10^9+7) 取模。

    解:

    可以看出,除了 (f(2)=3) 外,对于 (p e 2) ,均有(f(p)=p-1)

    同例题1,此时设 (g1(x)=x)(g2(x)=1)

    则此计算到结果时,除了 (f(1)) 没有计算外,还将 (f(2)) 的结果计算成 (1) ,所以最后答案还得加上 (3)

    #include<bits/stdc++.h>
    #define mem(a,b) memset(a,b,sizeof(a))
    using namespace std;
    typedef long long ll;
    const int maxn=1e6+50;
    const int mod=1e9+7;
    
    const int inv2=500000004;
    int pri[maxn],ncnt,T,m,id1[maxn],id2[maxn];
    bool isp[maxn];
    ll n,a[maxn],g1[maxn],g2[maxn],sum1[maxn];
    void init()
    {
        T=sqrt(n+0.5);
        for(int i=2;i<=T;i++){
            if(!isp[i])pri[++ncnt]=i,sum1[ncnt]=(sum1[ncnt-1]+i)%mod;
            for(int j=1;j<=ncnt&&i*pri[j]<=T;j++){
                isp[i*pri[j]]=true;
                if(i%pri[j]==0)break;
            }
        }
        for(ll l=1,r;l<=n;l=r+1){
            r=n/(n/l);
            a[++m]=n/l;
            if(a[m]<=T)id1[a[m]]=m;else id2[n/a[m]]=m;
            g1[m]=a[m]%mod;
            g1[m]=(g1[m]*(g1[m]+1)%mod*inv2%mod-1+mod)%mod;
            g2[m]=(a[m]-1+mod)%mod;
        }
        for(int i=1;i<=ncnt;i++){
            ll p2=1ll*pri[i]*pri[i],v;
            for(int j=1,id;j<=m&&p2<=a[j];j++){
                v=a[j]/pri[i];
                id=v<=T?id1[v]:id2[n/v];
                g1[j]=(g1[j]-1ll*pri[i]*(g1[id]-sum1[i-1]+mod)%mod+mod)%mod;
                g2[j]=(g2[j]-(g2[id]-(i-1)+mod)%mod+mod)%mod;
            }
        }
    }
    ll S(ll x,int y)
    {
        if(pri[y]>=x)return 0;
        int id=x<=T?id1[x]:id2[n/x];
        ll sum=(g1[id]-g2[id]-(sum1[y]-y+mod)%mod+mod+mod)%mod;
        for(int k=y+1;k<=ncnt&&1ll*pri[k]*pri[k]<=x;k++){
            ll pe=pri[k];
            for(int e=1;pe<=x;pe*=pri[k],e++)
                sum=(sum+1ll*(pri[k]^e)%mod*(S(x/pe,k)+(e!=1))%mod)%mod;
        }
        return sum;
    }
    int main()
    {
        scanf("%lld",&n);
        if(n==1){
            puts("1");return 0;
        }
        init();
        printf("%lld
    ",(S(n,0)+3)%mod);
    }
    



    例题3:2020ccpc网络赛 Graph Theory Class

    题意:

    (T) 个测试用例,(Tle 50) 。 对于每个测试用例给定的 (n,k) ,输出 $frac{n(n+3)}2+f(n+1)-4 $ 在模 (k) 下的结果。

    其中 $f(x)=sum_{1le ple x}p $。

    解:

    此时设 (g(x)=x) ,只需计算 (G(n)) 即为答案。

    #include<bits/stdc++.h>
    #define mem(a,b) memset(a,b,sizeof(a))
    using namespace std;
    typedef long long ll;
    const int maxn=1e6+50;
    
    int mod,inv2;
    ll kpow(ll a,ll b){
        ll ans=1;
        while(b){
            if(b&1)ans=ans*a%mod;
            a=a*a%mod;
            b>>=1;
        }
        return ans;
    }
    int pri[maxn],ncnt,T,m,id1[maxn],id2[maxn];
    bool isp[maxn];
    ll n,a[maxn],g1[maxn],sum1[maxn];
    void init()
    {
        m=ncnt=0;
        T=sqrt(n+0.5);
        for(int i=2;i<=T;i++){
            if(!isp[i])pri[++ncnt]=i,sum1[ncnt]=(sum1[ncnt-1]+i)%mod;
            for(int j=1;j<=ncnt&&i*pri[j]<=T;j++){
                isp[i*pri[j]]=true;
                if(i%pri[j]==0)break;
            }
        }
        for(ll l=1,r;l<=n;l=r+1){
            r=n/(n/l);
            a[++m]=n/l;
            if(a[m]<=T)id1[a[m]]=m;else id2[n/a[m]]=m;
            g1[m]=a[m]%mod;
            g1[m]=(g1[m]*(g1[m]+1)%mod*inv2%mod-1+mod)%mod;
        }
        for(int i=1;i<=ncnt;i++){
            ll p2=1ll*pri[i]*pri[i],v;
            for(int j=1,id;j<=m&&p2<=a[j];j++){
                v=a[j]/pri[i];
                id=v<=T?id1[v]:id2[n/v];
                g1[j]=(g1[j]-1ll*pri[i]*(g1[id]-sum1[i-1]+mod)%mod+mod)%mod;
            }
        }
    }
    
    int main()
    {
        int T;
        scanf("%d",&T);
        while(T--)
        {
            scanf("%lld%d",&n,&mod);
            inv2=kpow(2,mod-2);
            n++;init();n--;n%=mod;
            ll x=g1[id2[1]];
            printf("%lld
    ",(x-4+mod+n%mod*(n+3)%mod*inv2%mod)%mod);
        }
    }
    
  • 相关阅读:
    sql-trace-10046-trcsess-and-tkprof
    教你深入理解软件包的配置、编译与安装过程
    Java RESTful 框架的性能比较
    gcc、arm-Linux-gcc和arm-elf-gcc的组成及区别
    Linux线上系统程序debug思路及方法
    使用systemtap调试Linux内核 :www.lenky.info
    SystemTap使用技巧 1
    gvfs
    Systemtap examples, Network
    .NET 大型信息化建设标准基础数据管理平台
  • 原文地址:https://www.cnblogs.com/kkkek/p/13796551.html
Copyright © 2011-2022 走看看