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

    min_25筛是一个能快速求解积性函数前缀和的东西。

    要保证 $f(p)(p ext{是质数})$ 是个关于 $p$ 的多项式(次数也不要太高),并且 $f(p^k)$ 能快速计算。

    以下以洛谷的模板为例:($f(p^k)=(p^k)^2-p^k(p ext{是质数})$,求前 $N(Nle 10^{10})$ 项的和 $mod 10^9+7$ 的值)

    也就是 $sumlimits^N_{i=1}f(i)$。

    我们考虑分成质数和合数计算,同时枚举合数的最小质因子及其次数(因为是积性函数,所以直接提出来)。

    答案为 $sumlimits^N_{p=1}f(p)[pin prime]+sum f(p^e)[pin prime,p^ele N](sumlimits_{i=1}^{lfloorfrac{N}{p^e} floor}f(i)[i ext{的最小质因子}>p])$


    首先考虑质数部分:

    我们对每个次数不同的项分别计算。假设我们正在计算 $k$ 次方项。

    我们不知道为什么但就是设 $g(n,j)$ 为 $sumlimits^n_{i=1}i^k[i ext{是质数或者}i ext{的最小质因子}>p_j]$。

    其中 $p_j$ 是第 $j$ 小的质数。

    那么有初始状态 $g(n,0)=sumlimits^n_{i=2}i^k$(不把 $1$ 算进去)。

    转移,发现不满足 $g(n,j)$ 的限制且满足 $g(n,j-1)$ 的限制的数就是最小质因子恰为 $p_j$ 的合数。

    首先发现 $p_j^2>n$ 时这样的数不存在,此时 $g(n,j)=g(n,j-1)$。

    否则由于 $k$ 次方是完全积性函数,所以可以把 $p_j^k$ 提出来,那么这些数的 $k$ 次方和就是 $p_j^k(g(lfloordfrac{n}{p_j} floor,j)-g(p_{j-1},j-1))$。

    为什么前面的 $g$ 第二维是 $j$ 呢?因为我们只提出了一个 $p_j$,剩下可能还有。

    为什么要减掉后面的 $g$ 呢?因为后面是 $<p_j$ 的质数的 $k$ 次方和,如果把这一部分算上,那么就相当于整个答案多算了一些最小质因子小于 $p_j$ 的数。所以要减掉。(其实这一块就是质数 $k$ 次方的前缀和,可以预处理。以后记作 $s_k(j-1)$)

    最后用 $g(n,j-1)$ 减掉不满足 $g(n,j)$ 的限制且满足 $g(n,j-1)$ 的限制的数即可。

    方程出来了:

    $$g(n,j)=egin{cases}g(n,j-1)&p_j^2>n\g(n,j-1)-p_j^k(g(lfloordfrac{n}{p_j} floor,j)-s_k(j-1))&p_j^2le nend{cases}$$

    这个可以上滚动数组滚掉 $j$,同时对于 $p_j^2>n$ 的不用动就行了,时间复杂度也降到了……等等,这玩意复杂度明明比暴力还高好吗!

    但是我们可以发现,一个 $g(n,j)$ 可以由若干个 $g(lfloordfrac{n}{x} floor,k)$ 表示。

    这启示我们把所有 $lfloordfrac{N}{n} floor$ 相同的 $n$ 归为一类一起处理。具体实现,把 $g$ 的第一维下标变为 $lfloordfrac{N}{n} floor$。然而这样还是可能太大,于是可以对 $lesqrt{n}$ 和 $>sqrt{n}$ 的开两种编号。

    这一部分时间复杂度约为 $O(frac{N^{3/4}}{log N})$。

    怎么算的?我们发现类似埃氏筛,一个数 $x$ 被更新过 $frac{sqrt{x}}{logsqrt{x}}$ 次。

    那么复杂度 $O(sumlimits_{i=1}^{sqrt{N}}frac{sqrt{i}}{logsqrt{i}}+sumlimits_{i=1}^{sqrt{N}}frac{sqrt{frac{N}{i}}}{logsqrt{i}})$。经过积分爆算,发现略小于 $O(frac{N^{3/4}}{log N})$。

    先放一下模板的代码实现:

    inline int id(ll x){return x<=sq?id1[x]:id2[n/x];}
    void init(){
        sq=sqrt(n);
        FOR(i,2,sq){
            if(!vis[i]){
                pri[++pl]=i;
                s1[pl]=(s1[pl-1]+i)%mod;
                s2[pl]=(s2[pl-1]+(ll)i*i)%mod;    //质数一次项和二次项前缀和 
            }
            FOR(j,1,pl){
                if(i*pri[j]>sq) break;
                vis[i*pri[j]]=true;
                if(i%pri[j]==0) break;
            }
        }
        for(ll l=1,r;l<=n;l=r+1){
            r=n/(n/l);
            w[++tot]=n/l;
            if(n/l<=sq) id1[n/l]=tot;
            else id2[n/(n/l)]=tot;    //设置编号 
            int x=w[tot]%mod;
            g1[tot]=(1ll*x*(x+1)%mod*inv2%mod-1+mod)%mod; 
            g2[tot]=(1ll*x*(x+1)%mod*(2*x+1)%mod*inv6%mod-1+mod)%mod;    //2到x的和 
        }
    }
    void calc_g(){
        FOR(i,1,pl){    //遍历质数(上文中的j) 
            FOR(j,1,tot){    //遍历上文中的n 
                if((ll)pri[i]*pri[i]>w[j]) break;    //只管p_j^2<=n的 
                g1[j]=(g1[j]-1ll*pri[i]*(g1[id(w[j]/pri[i])]-s1[i-1]+mod)%mod+mod)%mod;
                g2[j]=(g2[j]-1ll*pri[i]*pri[i]%mod*(g2[id(w[j]/pri[i])]-s2[i-1]+mod)%mod+mod)%mod;
                //上面这两行就死磕吧……记得g1,g2和w下标的含义就好 
            }
        }
    }
    View Code

    再看整个式子:

    我们吸取了经验知道要绞尽脑汁地设 $S(n,j)$ 表示 $sumlimits_{i=2}^n f(i)[i ext{的最小质因子}>p_j]$。注意这里是 $f$ 而不是 $k$ 次方,下界是 $2$。

    那么要求即为 $S(N,0)+f(1)$。边界 $S(n,j)=0(p_j>n)$。

    转移,对答案有贡献的质数必定大于 $p_j$,那么贡献就是 $g(n,pl)-sum s_k(j)$。($pl$ 是 $lesqrt{N}$ 的质数个数)

    按照套路再枚举合数的最小质因子和次数,就是 $sum f(p^e)[pin prime,p>p_j,p^ele n](S(lfloordfrac{n}{p^e} floor)+[e e 1])$。

    (为何后面要有 $[e e 1]$ 呢?因为 $e=1$ 时,这个值在前面算质数时已经算过,而如果 $e e 1$,那么前面算质数时没算过,而因为 $S$ 下界为 $2$ 后面也不会算到,所以要加上 $f(p^e)$)。

    即使不记忆化,这一部分复杂度也是 $O(frac{N^{3/4}}{log N})$。(这个复杂度我是真不会分析了)

    贴上整题代码:

    #include<bits/stdc++.h>
    using namespace std;
    typedef long long ll;
    const int maxn=1000100,mod=1000000007,inv2=500000004,inv6=166666668;
    #define FOR(i,a,b) for(int i=(a);i<=(b);i++)
    #define ROF(i,a,b) for(int i=(a);i>=(b);i--)
    #define MEM(x,v) memset(x,v,sizeof(x))
    int sq,pri[maxn],pl,tot,g1[maxn],g2[maxn],id1[maxn],id2[maxn],s1[maxn],s2[maxn];
    bool vis[maxn];
    ll n,w[maxn];
    inline int add(int a,int b){return a+b<mod?a+b:a+b-mod;}
    inline int sub(int a,int b){return a<b?a-b+mod:a-b;}
    inline int mul(int a,int b){return (ll)a*b%mod;}
    inline int id(ll x){return x<=sq?id1[x]:id2[n/x];}
    void init(){
        sq=sqrt(n);
        FOR(i,2,sq){
            if(!vis[i]) pri[++pl]=i,s1[pl]=add(s1[pl-1],i),s2[pl]=add(s2[pl-1],mul(i,i));
            FOR(j,1,pl){
                if(i*pri[j]>sq) break;
                vis[i*pri[j]]=true;
                if(i%pri[j]==0) break;
            }
        }
        for(ll l=1,r;l<=n;l=r+1){
            r=n/(n/l);
            w[++tot]=n/l;
            if(n/l<=sq) id1[n/l]=tot;
            else id2[n/(n/l)]=tot;
            g1[tot]=sub(mul(mul(w[tot]%mod,(w[tot]+1)%mod),inv2),1);
            g2[tot]=sub(mul(mul(mul(w[tot]%mod,(w[tot]+1)%mod),(2*w[tot]+1)%mod),inv6),1);
        }
    }
    void calc_g(){
        FOR(i,1,pl){
            FOR(j,1,tot){
                if((ll)pri[i]*pri[i]>w[j]) break;
                g1[j]=sub(g1[j],mul(pri[i],sub(g1[id(w[j]/pri[i])],s1[i-1])));
                g2[j]=sub(g2[j],mul(mul(pri[i],pri[i]),sub(g2[id(w[j]/pri[i])],s2[i-1])));
            }
        }
    }
    int solve(ll nn,int x){    //S(nn,x)
        if(pri[x]>=nn) return 0;
        int ans=sub(sub(g2[id(nn)],s2[x]),sub(g1[id(nn)],s1[x]));    //题目中f为二次项减一次项 
        FOR(i,x+1,pl){    //p
            if((ll)pri[i]*pri[i]>nn) break;
            ll pro=1;
            FOR(j,1,50){    //e
                pro*=pri[i];    //p^e
                if(pro>nn) break;
                ans=add(ans,mul(mul(pro%mod,(pro-1)%mod),add(solve(nn/pro,i),j!=1)));
            }
        }
        return ans;
    }
    int main(){
        scanf("%lld",&n);
        init();
        calc_g();
        printf("%d
    ",add(solve(n,0),1));
    }
    View Code

    一些题目:

    LOJ6235 区间素数个数 题解

    LOJ6053 简单的函数 题解

  • 相关阅读:
    第07组 Alpha冲刺(1/6)
    第07组 团队Git现场编程实战
    第07组 团队项目-需求分析报告
    团队项目-选题报告
    第二次结对编程作业
    0012---求滑动距离
    0011---绝对值函数
    0010---温度转换
    0009---乘法问题
    0008---三位数倒序问题
  • 原文地址:https://www.cnblogs.com/1000Suns/p/10799229.html
Copyright © 2011-2022 走看看