Min_25筛学习Tip+链接

zhouyuheng2003

已于 2022-09-09 15:16:31 修改

阅读量1.9k

点赞数 2

分类专栏： OI 数论 Min_25筛文章标签：学习机器学习算法

于 2018-12-07 13:25:45 首次发布

本文链接：https://blog.csdn.net/zhouyuheng2003/article/details/84862000

版权

OI 同时被 3 个专栏收录

104 篇文章 2 订阅

订阅专栏

数论

14 篇文章 0 订阅

订阅专栏

Min_25筛

2 篇文章 0 订阅

订阅专栏

2022/9/9update：修缮及链接的更新

前言

机房里差不多都会Min_25筛了，我也赶紧补一波坑。
参考：txc的Min_25筛学习笔记，和yx的Min_25筛学习笔记。
~~前面两位dalao的都是过程详细的笔记，这里只记录我的一些思考，充其量就是一个小记。~~ 同时，2018的论文集讲的也很不错，内容规范，可以参考。

什么是Min_25筛

实际上，存在一种实现难度更低，同时在较小的数据范围下表现更好的方法，它和洲阁筛的思想类似，但理解起来更为容易
——朱震霆, 《一些特殊的数论函数求和问题》, 国家集训队2018论文集。

~~能代替洲阁筛，所以赶紧来学一波操作。~~
在不考虑常数，只考虑复杂度的情况下，Min_25筛是比杜教筛快的，Min_25的复杂度是 $\Theta(\frac{n^{\frac34}}{log\ n })$ 的，而杜教筛的复杂度是 $\Theta(n^{\frac23})$ ，使用Mathematica计算可以发现，当 $n$ 的值满足 $2.08977<n<3.00405*10^{22}$ 的时候，Min_25筛的复杂度是更优的。左端点…不用管，右端点的话，如果 $n$ 到了这个级别，按这两个复杂度已经做不了了(~~别想着打表，跑出这样的数据可能要跑100+天，再说空间复杂度也受不了~~)，所以Min_25筛的效率还是非常优越的。
Min_25筛空间复杂度为 $\Theta(\sqrt n)$ 。
再考虑适用范围：
杜教筛限制重重，需要更麻烦的推式子，而看Min_25适用的范围的介绍：

如果一个积性函数 f(i)在 i是质数时是一个关于 i的低次多项式，并且在质数的幂处的能快速求，那么大概可以用Min_25筛来求
——txc

相比之下，Min_25筛无疑是适用范围较广的
(~~yx:现在又找到杜教筛能做的而Min_25筛不能做的题吗~~，如果有欢迎评论)

知识储备

数论相关的基础知识（如果不太了解可以参考我的博客数论学习）

正文

想看推导请移步别人的链接，就不献丑了。
等博主学完Min_25筛后，发现别人写的比较清楚，我最多在后面给出一些给自己看的补充，就能完全理解Min_25筛了。
自己再写一遍没有什么意义。

学习指南（生存指南）

step1 打开zzt的论文
作为一名OIer你想必应该知道在哪里有论文。
zzt的论述还是很严密的，在思路上能给你很大的启发。
step2 打开txc的Min_25筛学习笔记
这一篇博客的思路是和zzq一样的，但是描述更加通俗，并且步骤更加简洁（没有冗长的式子）看完你应该能懂Min_25的更具体的原理与实现
注意：txc的博客中的一些定义回合zzt的论文中的定义不同，请仔细看，并且，对于第二部分，zzt的做法对应txc的第二种实现方式（似乎第二种方式常数小，但受到一些局限）
step3 看我写的内容
我是按这个步骤做下来的，由于我水平有限，所以看懂了式子后不知道怎么实现+不清楚复杂度，而看代码是没那么容易就清楚做法的，所以，在这里我会给一些Tip：
Tip 1：
在txc的博客中，他所说的“第一部分”的最后所说的记录 $\sqrt n$ 个值指的是记录当前已经求好的 $g (*, b - 1)$ ，现在要求的是 $g (*, b)$
容易发现，我们这里的 $b$ 的最大值是 $2$ 至 $\sqrt n$ 的质数数量，所以总的 $b$ 的取值是 $\Theta(\frac {\sqrt n}{ln\ {\sqrt n}})$ 的，这个值是可以接受的
考虑每一次转移，txc的博客里写到了 $a<prime^2$ 的都是不用转移的
剩下的就从最后一个值向前枚举更新，由于转移只会用到前面的值，所以只用开一个数组就够了，由于 $a$ 的取值是 $\left\lfloor\frac nm \right\rfloor$ 形式的，所以只用 $\Theta(\sqrt n)$ 种，空间是 $\Theta(\sqrt n)$ 的
Tip 2：
在写代码的时候有很多可以优化常数的地方，这个时候就要看仔细了，只有获取合适的写法才能让算法既好写又跑的快，并且对于Min_25筛，要注意哪些地方要用long long哪些地方不用，这是很容易错的（~~我自己就WA了好多发~~）
step4 自己写一遍Min_25筛，并做相关的习题
这里采用一道例题，BZOJ3944 Sum
这道题其实也是杜教筛模板题，这里用Min_25筛的板子也不错
贴出代码

#include<cstdio>
#include<cctype>
#include<cstring>
#include<algorithm>
#include<cmath>
namespace fast_IO
{
    const int IN_LEN=10000000,OUT_LEN=10000000;
    char ibuf[IN_LEN],obuf[OUT_LEN],*ih=ibuf+IN_LEN,*oh=obuf,*lastin=ibuf+IN_LEN,*lastout=obuf+OUT_LEN-1;
    inline char getchar_(){return (ih==lastin)&&(lastin=(ih=ibuf)+fread(ibuf,1,IN_LEN,stdin),ih==lastin)?EOF:*ih++;}
    inline void putchar_(const char x){if(oh==lastout)fwrite(obuf,1,oh-obuf,stdout),oh=obuf;*oh++=x;}
    inline void flush(){fwrite(obuf,1,oh-obuf,stdout);}
}
using namespace fast_IO;
#define getchar() getchar_()
#define putchar(x) putchar_((x))
#define rg register
typedef long long LL;
template <typename T> inline T max(const T a,const T b){return a>b?a:b;}
template <typename T> inline T min(const T a,const T b){return a<b?a:b;}
template <typename T> inline void mind(T&a,const T b){a=a<b?a:b;}
template <typename T> inline void maxd(T&a,const T b){a=a>b?a:b;}
template <typename T> inline T abs(const T a){return a>0?a:-a;}
template <typename T> inline void swap(T&a,T&b){T c=a;a=b;b=c;}
template <typename T> inline T gcd(const T a,const T b){if(!b)return a;return gcd(b,a%b);}
template <typename T> inline T lcm(const T a,const T b){return a/gcd(a,b)*b;}
template <typename T> inline T square(const T x){return x*x;};
template <typename T> inline void read(T&x)
{
    char cu=getchar();x=0;bool fla=0;
    while(!isdigit(cu)){if(cu=='-')fla=1;cu=getchar();}
    while(isdigit(cu))x=x*10+cu-'0',cu=getchar();
    if(fla)x=-x;
}
template <typename T> inline void printe(const T x)
{
    if(x>=10)printe(x/10);
    putchar(x%10+'0');
}
template <typename T> inline void print(const T x)
{
    if(x<0)putchar('-'),printe(-x);
    else printe(x);
}
LL sum1[700001];
int T,n,part,tot,qz[700001],prime[700001],sq[700001],primesize,sum0[700001],id[700001];
inline int ck(const int x){return x<=part?x:tot-n/x+1;}
int calc_mu(const int a,const int b)
{
    if(a<prime[b])return 0;
    int ans=sum0[ck(a)]+b-1;
    for(rg int i=b;i<=primesize&&sq[i]<=a;i++)ans-=calc_mu(a/prime[i],i+1);
    return ans;
}
LL calc_phi(const int a,const int b)
{
    if(a<prime[b])return 0;
    LL ans=sum1[ck(a)]-(qz[b-1]-(b-1));
    for(rg int i=b;i<=primesize&&sq[i]<=a;i++)
    {
        ans+=(calc_phi(a/prime[i], i+1)+prime[i])*(prime[i]-1);
        for(rg int x=prime[i]*prime[i],mine=x-prime[i];(LL)x*prime[i]<=a;x*=prime[i],mine*=prime[i])ans+=(calc_phi(a/x, i+1)+prime[i])*mine;
    }
    return ans;
}
int main()
{
    read(T);
    while(T--)
    {
        read(n);
        if(!n){print(0),putchar(' '),print(0),putchar('\n');continue;}
        part=sqrt(n);
        tot=primesize=0;
        for(rg int i=1;i<=part;i++)id[++tot]=i,sum0[tot]=i-1,sum1[tot]=((((LL)i+1)*i)>>1)-1;
        id[++tot]=n/part;if(id[tot]!=part)sum0[tot]=id[tot]-1,sum1[tot]=((LL)(((LL)id[tot]+1)*id[tot])>>1)-1;else tot--;
        for(rg int i=part-1;i>=1;i--)id[++tot]=n/i,sum0[tot]=id[tot]-1,sum1[tot]=((((LL)id[tot]+1)*id[tot])>>1)-1;
        for(rg int i=2;i<=part;i++)
            if(sum0[i]!=sum0[i-1])
            { 
                const int limit=(LL)i*i;
                for(rg int j=tot;id[j]>=limit;j--)
                {
                    const int t=ck(id[j]/i);
                    sum1[j]-=(sum1[t]-qz[primesize])*i;
                    sum0[j]-=sum0[t]-primesize;
                }
                prime[++primesize]=i,sq[primesize]=i*i,qz[primesize]=qz[primesize-1]+i;
            }
        for(rg int i=1;i<=tot;i++)sum1[i]-=sum0[i],sum0[i]=-sum0[i];//sum1为phi的前缀和，sum0为mu的前缀和
        print(calc_phi(n,1)+1),putchar(' '),print(calc_mu(n,1)+1),putchar('\n');
    }
    return flush(),0;
}