zoukankan      html  css  js  c++  java
  • min_25筛 学习笔记

    min_25筛 学习笔记

    简介

    min_25筛是一种筛法,由min_25提出,能够 (O(dfrac{n^{frac{3}{4}}}{ln n})) 的求积性函数 (f) 的前缀和。而 (f) 需要满足如下几个条件:

    (设 (p) 为任意质数)

    • (f(p)) 可以写成低次多项式的形式。设这个多项式为 (g(x))
    • (f(p^k)) 很方便计算

    约定

    要求前缀和的积性函数为 (f),在质数位置可以写成多项式 (g)

    设第 (i) 个质数为 (p_i)。设 ([1,n]) 中质数总共有 (|p|) 个。

    (n) 的“基本和组”定义为集合 (B(n)={x|x=n/k, kin[1,n]})。我们知道 (|B(n)|) 大约是 (2sqrt{n}),根据整除分块。

    求法

    part-1 质数位置

    首先我们求出 ([1,n]) 中质数位置的 (f) 的和,也就是 (g) 的和。

    (g) 是一个低次的多项式,因此我们可以把它的每一次单独拿出来做。那我们相当于要求 ([1,n]) 中的质数的 (k) 次幂和。

    考虑埃氏筛的过程。每次我们选一个质数 (p),并且把它的倍数筛去。假设我们只看 “有效筛” 的位置:即,如果这个位置已经被筛掉,就不去动它。如果 (x) 被“有效筛” 了,那 (x) 最小的质因子是 (p),且 (x) 为合数。那么 (xge p^2)

    也就是说,这个过程中,我们只需要 (le sqrt{x}) 的那些质数。这些质数是可以预先筛出来的。

    我们记下每次筛完之后的结果。换句话说,记 (G(n,i)) 表示,([1,n]) 中,只留下素数与最小质因子 (> p_i) 的合数,这些数的 (g) 的和。

    考虑 (G(n,i-1) o G(n,i)),相当于要把 (i) 能 “有效筛” 的那些位置的贡献去掉。

    (x) 能被“有效筛”,(xle n)。那 (x) 一定是 (p_i) 的倍数,且 (x/p_i) 的最小质因子 (ge p_i),也就是说,(> p_{i-1})

    注意到 (xle n),那 (x/p_ile n/p_i)。且 (x) 的最小质因子 (>p_{i-1})。满足条件的这些 (x) 大约就是 (G(n/p_i,i-1))。但是,这个 (G) 还把质数算上去了,我们要抠掉比 (p_i) 小的质数。 (注意到,如果等于 (p_i) 是没问题的,这相当于 (x=p_i^2),满足的)

    这相当于前缀若干个质数的 (g) 的和。可以在筛质数的时候顺便求个前缀和得到,记 (sump(i)) 表示前 (i) 个质数的 (g) 和。

    那么我们得到

    [Delta=g(p_i) imes left[G(n/p_i,i-1)-sump(i-1) ight] ]

    其中 (Delta) 表示变化量,即 (G(n,i-1)-G(n,i))。 也就是说我们令 (G(n,i)=G(n,i-1)-Delta) 即可。

    边界:(G(n,0)=sumlimits_{i=2}^{n} g(k)),把 (g) 的每一项拆开,它就可以直接计算了。

    备注:

    我们在实现的时候,对 (i) 这一维滚动,且只保留基本和组位置的值。

    由于后面只需要 (G(*,tot)) 位置的值,所以后面直接记

    [G(i)=sumlimits_{ple i,pin prime} g(p) ]

    可以用积分证明,这部分复杂度是 (O(dfrac{n^{frac{3}{4}}}{ln n})) 的。

    part-2 求答案

    我们已经知道了质数位置的 (g) 和。

    类似的考虑,设 (S(n,i)) 表示,最小质因子 (ge p_i) 的那些数的 (f) 的和。

    考虑一下边界:(S(n,|p|+1)=0) 。显然,因为 ([1,n]) 中的数,最小质因子都 (<p_{|p|})

    然后我们要求的答案是 (S(n,1)+f(1))。这也显然,因为 (1) 以上所有数的最小质因子都 (ge 2)。那 (S(n,1)) 就会覆盖到 (2) 以上所有数。

    注意到我们是知道一个大值,然后要求的是一个小值。那考虑从大往小推。即,假设大的已知,做小的。

    假设现在要做 (S(n,i))

    首先加上满足条件质数的答案。这玩意就是 ([1,n]) 中质数的 (g) 和,减去 (<p_i) 的质数的 (g) 和。即,(G(n)-sump(i-1))

    然后是满足条件合数的答案。如何钦点合数呢?合数就至少有两个质因数。枚举它的最小质因数是谁,然后枚举它多少次,然后枚举后面选了啥数。

    设枚举的最小质因数是 (j (jge i)),然后 (p_j)(e) 次方。那它后面可能还跟了一个 (p_j),也可能跟的是比 (p_j) 大的质因数。

    注意到 (f(p_j^{e+1})) 需要特判一下,其它的可以用 (S) 写出来。脑子转一下,得到它的贡献为

    [f(p_j^{e+1})+f(p_j^{e}) imes S(n/p_j^{e},i+1) ]

    为啥后面直接乘起来是对的,因为 (f) 是积性函数。

    需要满足:(p_j^{e+1}le n)。这样枚举一遍,把贡献加起来即可。

    最后 (S(n,1)+1) 就是答案力!!1

    板子题

    loj6053 简单的函数

    分析一下,发现它在质数的时候可以写成 (g(p)=p-1) 的形式,特判 (2)

    然后就是个多项式,用 min-25 筛就可以做了。

    代码中有注释,可以参考一下实现。

    #include <bits/stdc++.h>
    using namespace std;
    namespace Flandre_Scarlet
    {
        #define N   200005
        #define int long long
        #define mod 1000000007
        #define iv2  500000004
        #define F(i,l,r) for(int i=l;i<=r;++i)
        #define D(i,r,l) for(int i=r;i>=l;--i)
        #define Tra(i,u) for(int i=G.h[u],v=G.to(i);~i;i=G.nx(i),v=G.to(i)) if (i>=0)
        #define MEM(a,x) memset(a,x,sizeof(a))
        #define FK(a) MEM(a,0)
        #define sz(x) ((int)x.size())
        #define all(x) x.begin(),x.end()
        #define p_b push_back
        #define pii pair<int,int>
        #define fir first
        #define sec second
        int I() {char c=getchar(); int x=0; int f=1; while(c<'0' or c>'9') f=(c=='-')?-1:1,c=getchar(); while(c>='0' and c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar(); return ((f==1)?x:-x);}
        template <typename T> void Rd(T& arg){arg=I();}
        template <typename T,typename...Types> void Rd(T& arg,Types&...args){arg=I(); Rd(args...);}
        void RA(int *p,int n) {F(i,1,n) *p=I(),++p;}
        int pr[N]; bool notp[N];
        int sp[N];
        void Init() // 筛到sqrt(n)
        {
            int n=2e5;
            int&c=pr[0];
            notp[1]=1;
            F(i,2,n)
            {
                if (!notp[i])
                {
                    pr[++c]=i;
                }
                for(int j=1;j<=c and i*pr[j]<=n;++j)
                {
                    int u=pr[j];
                    notp[i*u]=1;
                    if (i%u==0) break;
                }
            }
            F(i,1,c) sp[i]=sp[i-1]+pr[i];
        }
        int n;
        void Input()
        {
            n=I();
        }
        
        int g1[N],g0[N]; // 1次方, 0次方的和
        int w[N],i1[N],i2[N];
        // w: 存储基本和组位置; i1,i2: 给这些位置编号
        // 若x<=sn,i1[x]保存x的编号; 否则, i2[n/x]保存x的编号
        int sn,tot;
        int s1(int x){x%=mod; return x*(x+1)%mod*iv2%mod;}
        int S(int x,int p) // [1,x], min-p>=pr[p], f sum
        {
            if (x<=1 or pr[p]>x) return 0;
            int k=(x<=sn)?(i1[x]):(i2[n/x]);
            int ans=(g1[k]-g0[k])-(sp[p-1]-(p-1)); // 质数部分的答案
            ans=(ans%mod+mod)%mod;
            if (p==1) ans+=2; // 特判2: 我们把f(2)当成了1, 实际上它是3
    
            for(int i=p;i<=pr[0] and pr[i]<=x/pr[i];++i)
            {
                int u=pr[i];
                int e=1,ue=u;
                // ue为u的e次方
                while(ue*u<=x) // u的e+1次方<=x
                {
                    ans+=((u^(e+1))+S(x/ue,i+1)*(u^e)%mod)%mod;
                    // u^e为 f(u的e次方), u^(e+1)同理
                    ans%=mod;
                    ue*=u; ++e;
                }
            }
            return (ans%mod+mod)%mod;
        }
        void Sakuya()
        {
            sn=sqrt(n),tot=0;
            for(int l=1,r;l<=n;l=r+1)
            {
                r=(n/(n/l));
                w[++tot]=n/l;
                if (n/l<=sn) i1[n/l]=tot;
                else i2[n/(n/l)]=tot;
    
                g0[tot]=((n/l)-1)%mod;
                g1[tot]=(s1(n/l)-1)%mod;
                // G(x)初始为g(2)+g(3)...+g(x)
                // 这里拆成两个, 因此 g0=x-1, g1=2+3...+(x-1)
            }
            F(j,1,pr[0])
            {
                int u=pr[j];
                for(int i=1;i<=tot and u*u<=w[i];++i)
                {
                    int t=w[i]/u,k;
                    if (t<=sn) k=i1[t];
                    else k=i2[n/t];
                    g0[i]-=g0[k]-(j-1);
                    g1[i]-=(g1[k]-sp[j-1])%mod*u%mod;
                    g0[i]=(g0[i]%mod+mod)%mod;
                    g1[i]=(g1[i]%mod+mod)%mod;
                }
            }
            // 计算G
            int ans=S(n,1)+1; // f(1)=1, 加上
            printf("%lld
    ",ans);
        }
        void IsMyWife()
        {
            Init();
            Input();
            Sakuya();
        }
    }
    #undef int //long long
    int main()
    {
        Flandre_Scarlet::IsMyWife();
        return 0;
    }
    
  • 相关阅读:
    新一代MQ apache pulsar的架构与核心概念
    Flutter使用fluwx实现微信分享
    BZOJ3622 已经没有什么好害怕的了 动态规划 容斥原理 组合数学
    NOIP2016提高组Day1T2 天天爱跑步 树链剖分 LCA 倍增 差分
    Codeforces 555C Case of Chocolate 其他
    NOIP2017提高组Day2T3 列队 洛谷P3960 线段树
    NOIP2017提高组Day2T2 宝藏 洛谷P3959 状压dp
    NOIP2017提高组Day1T3 逛公园 洛谷P3953 Tarjan 强连通缩点 SPFA 动态规划 最短路 拓扑序
    Codeforces 873F Forbidden Indices 字符串 SAM/(SA+单调栈)
    Codeforces 873E Awards For Contestants ST表
  • 原文地址:https://www.cnblogs.com/LightningUZ/p/15386158.html
Copyright © 2011-2022 走看看