menu Shadowice!
复健笔记min25筛
153 浏览 | 2020-09-24 | 阅读时间: 约 2 分钟 | 分类: ACM | 标签:

简介

min_25筛是一种用来求解积性函数前缀和的算法,分为两部分,第一部分复杂度为$O(\frac{n^{0.75}}{\log n})$,第二部分复杂度为$O(n^{1-\epsilon})$,也就是近线性复杂度,实际测试时可以通过$10^{10}-10^{12}$的数据

前置芝士

积性函数

基本的和式推演能力

记号与约定

$P_{j}$:第$j$大的质数

$ f(x) $ :待求前缀和的积性函数

$ g(x) $ :$f(x)$ 的取值在$x$为质数时等于另外一个函数$g(x)$

$mindiv(x)$:$x$的最小质因子

定义$pre(n)=\sum_{i=1}^{n}f(p_{i})$

定义$dp(n,j)=\sum_{p \leq n}g(p)$,满足$p$是质数或者$mindiv(p)>P_{j}$

定义$sum(n,j)=\sum_{2 \leq i \leq n,mindiv(i) \geq P_{j}}f(i)$

递归式

$ dp(n,j)=dp(n,j-1)-g(P_{j})(sp(\lfloor \frac{n}{P_j}\rfloor,j-1)-pre(j-1))$

$dp(n,0)=\sum_{i=2}^{n}g(i)$

$sum(n,j)=sp(n,\min_{P_j >\lceil \sqrt{n} \rceil}j)-pre(j-1)+\sum_{k \geq j}\sum_{e}f(P_{k}^{e})(sum( \lfloor \frac{n}{p_{k}^{e}} \rfloor,k+1)+[e>1])$

$sum(n,j)=0$,当$n < P_{j}$时成立。

$sum(n,j)=sp(n,\min_{P_j >\lceil \sqrt{n} \rceil}j)-pre(j-1),当P_{j}^2 >n 时成立$

算法实现

计算$sum$函数

为了计算$f$的前缀和,我们定义了$sum$函数,不难发现$\sum_{i=1}^{n}f(i)=sum(n,1)+1$,因此计算$sum(n,1)$即可完成任务

$sum $的递归式是这样推导的:前一部分是所有质数对sum的贡献,后一部分是所有合数对sum的贡献

我们知道当$P^2_{j}>n$时,不会存在最小质因子大于$P_j$的合数

那么此时的$dp(n,j)$就是$f(x)$在x取所有小于n的质数时的函数值之和,减掉比$P_j$小的质数值,剩下的就是符合要求的质数

所以$sum$递归式的前一部分就非常容易理解了

后一部分则是枚举合数的最小质因子$P_{k}$以及这个质因子的次幂$e$,然后做一个递归

这里着重讲一下$[e>1]$这一项的意义,因为$sum$的定义不包含1,因此直接递归我们会落下形如$P_{k}^{e},e>1$这样的合数,而当$e=1$时后面的sum就没必要再加$1$,因为质数的贡献已经在前一部分计算过了

计算$dp$函数

容易发现要用到的$n$值只有$\sqrt{n}$种,要用到的$j$值只有$\frac{\sqrt{n}}{\log n}$种,不过直接实现还是很慢

做一个剪枝,首先我们使用滚动数组,滚掉第二维j,完成这个dp

当$P_{j}^{2}>n$时,$dp(n,j),dp(n,j-1)$的值都是质数函数值的和,所以我们可以不用转移这些位置的dp值

此时复杂度被优化到了$O(\frac{n^{0.75}}{\log n})$可以满足我们的需要

适用范围

$g(x)$的前缀和能快速计算

$f(P_{k}^{e})$能快速计算

代码

下面给出了一个用min_25筛计算$\phi$函数前缀和代码示例,这份代码经过了良好的封装,只要改一下$sg,g,pre$函数就可以套用到别的题目当中

#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
namespace min_25
{
    typedef long long ll;
    const int N=1e5+1000;
    int zhi[N];int book[N];int ct;int ctcur;
    ll n;ll lim;
    ll q[N<<2];int hd;
    ll pre1[N];ll pre2[N];
    ll pre[N];
    ll g1(ll n)// 质数 
    { 
        return n;
    }
    ll sg1(ll n)// 前缀和 
    {
        return n*(n+1)/2;
    }
    ll g2(ll n){return 1;}
    ll sg2(ll n){return n;}
    ll fpw(ll pw,int k,int e)// pke 
    {
        return (ll)(zhi[k]-1)*pw/zhi[k];
    }
    void euler()
    {
        for(int i=2;i<=N-10;i++)
        {
            if(book[i]==0)zhi[++ct]=i;
            for(int j=1;j<=ct&&i*zhi[j]<=N-10;j++)
                {    
                    book[i*zhi[j]]=true;
                    if(i%zhi[j]==0)break;
                }
        }
        for(int i=1;i<=N-10;i++)
            pre1[i]=pre1[i-1]+g1(zhi[i]);
        for(int i=1;i<=N-10;i++)
            pre2[i]=pre2[i-1]+g2(zhi[i]);
        for(int i=1;i<=N-10;i++)
            pre[i]=pre1[i]-pre2[i];
    }
    void init(ll rn)
    {
        n=rn;lim=sqrt(n);
        ctcur=ct;
        while(zhi[ctcur-1]>lim)ctcur--;
        hd=0;
        for(ll i=1;i<=n;i=n/(n/i)+1)
            q[++hd]=n/i;    
    }
    struct shrink_array
    {
        ll a1[N];ll a2[N];
        ll& operator [](int x)
        {
            return (x<=lim)?a1[x]:a2[n/x];
        }
    }dp1,dp2,dp;
    void dp1_go()
    {
        for(int i=1;i<=hd;i++)
            dp1[q[i]]=sg1(q[i])-g1(1);
        for(int j=1;j<=ctcur&&(ll)zhi[j]*zhi[j]<=n;j++)
            for(int i=1;q[i]>=(ll)zhi[j]*zhi[j];i++)
                dp1[q[i]]-=g1(zhi[j])*(dp1[q[i]/zhi[j]]-pre1[j-1]);
    }
    void dp2_go()
    {
        for(int i=1;i<=hd;i++)
            dp2[q[i]]=sg2(q[i])-g2(1);
        for(int j=1;j<=ctcur&&(ll)zhi[j]*zhi[j]<=n;j++)
            for(int i=1;q[i]>=(ll)zhi[j]*zhi[j];i++)
                dp2[q[i]]-=g2(zhi[j])*(dp2[q[i]/zhi[j]]-pre2[j-1]);
    }
    
    ll sum(ll n,int j)
    {
        if(n<zhi[j])return 0;
        ll ret=dp[n]-pre[j-1];
        if(j>ctcur)return (n<=(ll)zhi[ctcur])?0:ret;
        for(int k=j;(ll)zhi[k]*zhi[k]<=n&&k<=ctcur;k++)
            for(ll e=1,pw=zhi[k];pw<=n;e++,pw*=zhi[k])
                ret+=fpw(pw,k,e)*(sum(n/pw,k+1)+(e>1));
        return ret;
    }
    ll solve(ll rn)
    {
        init(rn);
        dp1_go();
        dp2_go();
        for(int i=1;i<=hd;i++)
            dp[q[i]]=dp1[q[i]]-dp2[q[i]];
        return sum(n,1)+1;
    }
    void clear()
    {
        for(int i=1;i<=hd;i++)
            dp1[q[i]]=dp2[q[i]]=dp[q[i]]=0;
    }
}
int main()
{
    min_25::euler();
    int t;
    scanf("%d",&t);
    for(int i=1;i<=t;i++)
    {
        ll n;
        scanf("%lld",&n);
        printf("%lld\n",min_25::solve(n));
        min_25::clear();    
    }
    return 0;
}
    
知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议

发表评论

email
web

全部评论 (共 2 条评论)

    2020-09-29 19:29
    第二部分貌似加个记忆化复杂度就可以有$O(\frac{n^\frac{3}{4}}{\log n})$,虽然年代久远都点记不清了...
      shadowice1984
      2020-09-30 20:54
      @Mr_Spade应该是不行的,第二部分记忆化了反而更慢而且复杂度也不会变优秀,虽然年代久远我的记忆也已经模糊了