杜教筛小结

首先的首先，%%%%%%杜教

……咳咳，言归正传，现在来总结一下杜教筛……

杜教筛大概是一种用于快速处理积性函数前缀和的有力工具，主要利用了狄利克雷卷积的知识

我们用高中均值不等式的复杂度分析可以得到其复杂度是$O(n^{\frac{2}{3}})$的

具体的讲解我就不写了……留几篇比较好的博客。

https://www.cnblogs.com/abclzr/p/6242020.html

这里总结一下我做过的题目吧……

1.bzoj3944

由于太水以至于没打……

可以当做板子练习题的题目。

2.bzoj4176

深刻的记住了一个式子……之前一直没有这种意识，有一些转换的确很显然但是就是没想过

\[d\mid ij\Leftrightarrow \frac{d}{(i,d)}\mid j\]

……自己打公式太麻烦了，放弃……

有两种解法，一种是类似结论的公式，一种是暴力反演

……其实怎么演都能演出来，关键就看第一步转化能不能出来

代码：

#include <cstdio>
#include <map>
#include <cstring>
using namespace std;
#define mod 1000000007
#define N 1000000
#define LL long long
int n,prime[N/10],tot,vis[N+10],mu[N+10];
inline void intn()
{
    register int i,j,tmp;
    for(mu[1]=1,i=2;i<=N;++i)
    {
        if(!vis[i])prime[++tot]=i,mu[i]=-1;
        for(j=1;j<=tot&&(tmp=i*prime[j])<=N;++j)
            {vis[tmp]=1;if(i%prime[j]==0){mu[tmp]=0;break;}mu[tmp]=-mu[i];}
    }
    for(i=1;i<=N;++i)mu[i]+=mu[i-1];
}
map<int,int> mmp;
inline int getmu(int x)
{
    if(x<=N)return mu[x];
    if(mmp.count(x))return mmp[x];
    register int i,last,ret=1;
    for(i=2;i<=x;i=last+1)
        last=x/(x/i),ret=(ret-((LL)(last-i+1)*getmu(x/i)%mod))%mod;
    return mmp[x]=ret;
}
inline int calc2(int x)
{
    register int i,last,ret=0;
    for(i=1;i<=x;i=last+1)
        last=x/(x/i),ret=(ret+((LL)(last-i+1)*(x/i)%mod))%mod;
    return (LL)ret*ret%mod;
}
int main()
{
    // freopen("Ark.in","r",stdin);
    register int i,last,ans=0;
    scanf("%d",&n);intn();
    for(i=1;i<=n;i=last+1)
        last=n/(n/i),ans=(ans+((LL)(getmu(last)-getmu(i-1))*calc2(n/i)%mod) )%mod;
    printf("%d\n",(ans+mod)%mod);
}

3.bzoj3512

然后这题我们发现n比较小，m却比较大

所以考虑与枚举n有关的算法……

设$S(n,m)=\sum _{i=1} ^{n} \varphi (i*n)$

那么我们知道$\varphi$是个积性函数……所以上来质因数分解一发

设$n=\prod p_{i} ^{a_{i}}$

那么有$S(n,m)=p_{i}^{a_{i}-1} * S(\frac{n}{p_{i}^{a_{i}-1}},m)$

然后我们把p都干下去了之后，

考虑$\varphi(p*q)=\varphi(p)*\varphi(q),(p,q)==1$

那么我们有$S(n,m)=\varphi(p_{i}) *   \sum _{i=1} ^{m} \varphi(\frac{n}{p_{i}}*i)*[(i,p_{i})==1]+p_{i}* \sum _{i=1} ^{m} \varphi(n*\frac{i}{p_{i}})*[p_{i}|i]$

又因为$\varphi(p)==p-1$

所以我们可以凑个整……

$S(n,m)=\varphi(p_{i}) * \sum _{i=1} ^{m} \varphi(\frac{n}{p_{i}}*i)+\sum _{i=1} ^{ \left \lfloor \frac{m}{p_{i}} \right \rfloor } \varphi(n*i)$

$S(n,m)=S(\frac{n}{p_{i}},m)+S(n, \left \lfloor \frac{m}{p_{i}} \right \rfloor)$

这样我们就可以记忆化搞了

边界条件$n==1$就是$\varphi$前缀和，杜教筛之。

#include <cstdio>
#include <cstring>
#include <map>
using namespace std;
#define N 1000000
#define LL long long
#define mod 1000000007
int prime[N/10],tot,mind[N+10],phi[N+10],sum[N+10];
inline void intn()
{
    register int i,j,tmp;
    phi[1]=1;
    for(i=2;i<=N;++i)
    {
        if(!mind[i])prime[++tot]=i,mind[i]=i,phi[i]=i-1;
        for(j=1;j<=tot&&(tmp=i*prime[j])<=N;++j)
        {
            mind[tmp]=prime[j];
            if(i%prime[j]==0){phi[tmp]=phi[i]*prime[j];break;}
            phi[tmp]=phi[i]*phi[prime[j]];
        }
    }
    for(i=1;i<=N;++i)sum[i]=(sum[i-1]+phi[i])%mod;
}
inline int quick_mod(int di,int mi)
{
    int ans=1;
    for(di%=mod;mi;mi>>=1,di=(LL)di*di%mod)
        if(mi&1)ans=(LL)ans*di%mod;
    return ans;
}
inline int divide(int x)
{
    register int ret=1,i;
    for(i=1;(LL)prime[i]*prime[i]<=x;++i)
        if(x%prime[i]==0)
        {
            ret*=prime[i];
            while(x%prime[i]==0)x/=prime[i];
        }
    return ret*x;
}
#define N1 100010
map<int,int>mmp;
int mem[N1];
inline int getphi(int x)
{
    if(x<=N)return sum[x];
    else if(mmp.count(x))return mmp[x];
    int ret=((LL)x*(x+1)/2%mod);
    register int i,last;
    for(i=2;i<=x;i=last+1)
        last=x/(x/i),ret=( ret+mod-( (LL)(last-i+1)*getphi(x/i)%mod ) )%mod;
    return mmp[x]=ret;
}
inline int calc(int a,int b)
{
    if(!b)return 0;
    if(b==1)return phi[a];
    if(a==1)return getphi(b);
    return ((LL)calc(a/mind[a],b)*(mind[a]-1)%mod+calc(a,b/mind[a]))%mod;
}
int main()
{
    // freopen("Ark.in","r",stdin);
    register int n,m,ans=0,tmp,x,i;
    scanf("%d%d",&n,&m),intn();
    memset(mem,-1,sizeof(mem));
    for(i=1;i<=n;++i)
    {
        x=divide(i),tmp=i/x;
        if(mem[x]==-1)mem[x]=calc(x,m);
        ans=(ans+(LL)tmp*mem[x]%mod)%mod;
    }
    printf("%d\n",ans);
}

4.bzoj4916

 首先你发现$\sum \mu(i^{2})$是出来搞笑的，输出1即可

然后我们考虑一下这个$\sum _{i=1} ^{n} \varphi(i^{2})$

由于$\varphi(i)$是积性函数，所以我们考虑

$\varphi(i)=\prod \varphi(p_{i})*p_{i} ^{a_{i}-1}$

$\varphi(i^{2})=\prod \varphi(p_{i})*p_{i} ^{a_{i}-1+a_{i}}$

$\varphi(i^{2})=\prod \varphi(p_{i})*p_{i} ^{a_{i}-1}   *p_{i}^{a_{i}}$

所以……$\varphi(i^{2})=i*\varphi(i)$

设$f(i)=\varphi(i^{2})$

那么我们不是知道$\sum _{d|n} \varphi(d)  = n$ 嘛

那么转化一下 $\sum _{d|n} d*\varphi(d) * \frac {n}{d} = n * \sum _{d|n} \varphi(d) =n^{2}$

所以我们有积性函数$f(i)=i\varphi(i),g(i)=i,h(i)=i^{2}$

这样就可以杜教筛了，设大写字母为前缀和，过程略……

最后结果是$F(i)=H(i)-\sum _{j=2} ^{m} g(j)*F(\left \lfloor \frac{m}{j} \right \rfloor )$

其中，$G(i)=\frac {n*(n+1)} {2} , H(i)= \frac {n*(n+1)*(2n+1)} {6} $

代码：

#include <cstdio>
#include <map>
#include <cstring>
using namespace std;
#define N 1000000
#define mod 1000000007
int n,inv6,vis[N+10],prime[N/10],tot,phi[N+10],sum[N+10];
map<int,int> mmp;
#define LL long long
inline void intn()
{
    register int i,j,tmp;
    for(phi[1]=1,i=2;i<=N;++i)
    {
        if(!vis[i])prime[++tot]=i,phi[i]=i-1;
        for(j=1;j<=tot&&(tmp=i*prime[j])<=N;++j)
        {
            if(i%prime[j]==0)
                {vis[tmp]=1,phi[tmp]=phi[i]*prime[j];break;}
            vis[tmp]=1,phi[tmp]=phi[i]*phi[prime[j]];
        }
    }
    for(i=1;i<=N;++i)
        sum[i]=(sum[i-1]+(LL)i*phi[i]%mod)%mod;
}
inline int quick_mod(int di,int mi)
{
    int ans=1;
    for(;mi;mi>>=1,di=(LL)di*di%mod)
        if(mi&1)ans=(LL)ans*di%mod;
    return ans;
}
inline int getf(int x)
{
    if(x<=N)return sum[x];
    if(mmp.count(x))return mmp[x];
    register int i,last,ret=(LL)x*(x+1)%mod*(2*x+1)%mod*inv6%mod;
    for(i=2;i<=x;i=last+1)
        last=x/(x/i),
          ret=(   ret +mod -(LL)getf(x/i)*( (LL)(last-i+1)*(last+i)/2 %mod )%mod   )%mod; 
    return mmp[x]=ret;
}
int main()
{
    // freopen("Ark.in","r",stdin);
    inv6=quick_mod(6,mod-2);
    scanf("%d",&n),intn();
    printf("1\n%d\n",getf(n));
}

5.bzoj3930

我们先反演一波……

定义$f(i)$为$gcd(a_{1},a_{2},......,a_{n})==i$的方案数,$F(i)$为$i|gcd(a_{1},a_{2},......,a_{n})$的方案数

那么有$F(i)=\sum _{i|n} f(n)$

然后我们演他一下，可以得到：

$f(n)=\sum _{n|d} \mu(\frac{d}{n})*F(d)$

我们考虑一下……$F(d)=( \left \lfloor \frac{r}{d} \right \rfloor -  \left \lfloor \frac{l-1}{d} \right \rfloor )^{n}$

然后我们考虑枚举$i=\frac{d}{n}$

$f(n)=\sum_{i=1}^{  \left \lfloor \frac{r}{n} \right \rfloor }  \mu(i) * ( \left \lfloor \frac{r}{i*n} \right \rfloor -  \left \lfloor \frac{l-1}{i*n} \right \rfloor )^{n}$

然后我们对$\left \lfloor \frac{r}{i*n} \right \rfloor$和$\left \lfloor \frac{l-1}{i*n} \right \rfloor$除法分块

再杜教筛一个$\mu$的前缀和就行了。

代码：

#include <cstdio>
#include <map>
#include <cstring>
using namespace std;
#define mod 1000000007
#define LL long long
#define inf 0x3f3f3f3f
#define N 1000000
int prime[N],tot,mu[N+10],vis[N+10];
int n,d,l,r;
inline void intn()
{
    register int i,j,tmp;
    for(mu[1]=1,i=2;i<=N;++i)
    {
        if(!vis[i])prime[++tot]=i,mu[i]=-1;
        for(j=1;j<=tot&&(tmp=i*prime[j])<=N;++j)
        {
            vis[tmp]=1;
            if(i%prime[j]==0)break;
            mu[tmp]=-mu[i];
        }
    }
    for(i=2;i<=N;++i)mu[i]+=mu[i-1];
}
inline int min(int a,int b){return a<b?a:b;}
map<int,int> mmp;
inline LL getmu(int x)
{
    if(x<=N)return mu[x];
    if(mmp.count(x))return mmp[x];
    LL ret=1;
    for(register int i=2,last;i<=x;i=last+1)
        last=x/(x/i),ret=( ret-getmu(x/i)*(LL)(last-i+1)%mod )%mod;
    return mmp[x]=ret;
}
inline LL quick_mod(LL di,int mi)
{
    LL ans=1;
    for(;mi;mi>>=1,di=di*di%mod)
        if(mi&1)ans=ans*di%mod;
    return ans;
}
signed main()
{
    register int i,last;
    scanf("%d%d%d%d",&n,&d,&l,&r);
    r=r/d,l=(l-1)/d,intn();
    LL ans=0;
    for(i=1;i<=r;i=last+1)
    {
        last= (l/i) ? min( r/(r/i), l/(l/i) ) : r/(r/i) ,
        ans=(ans+ ( getmu(last)-getmu(i-1) )*quick_mod((r/i-l/i),n) )%mod;
    }
    printf("%lld\n",(ans%mod+mod)%mod);
}

6.bzoj4652

…………不得不说是个好题

被拿去做考试题了

结果只打了24分，而不是常见的84分暴力……原因待会再说

介绍一下考试历程

首先我观察了一下十进制下都以谁为分母是纯循环的……

发现2，5，10都不是

然后发现1无论何时是纯循环的

然后发现这是个完全积性函数！！！哇太棒了！

于是我就把这个函数命名为$f(i)$，当$\frac{1}{i}$纯循环时$f(i)$为1，否则为0

然后我就开始反演……

$\sum _{i=1} ^{n} \sum _{j=1} ^{m} f(j)[(i,j)==1]$

$\sum _{i=1} ^{n} \sum _{j=1} ^{m} f(j) \sum _{d|(i,j)}\mu(d)$

然后我们枚举$d$得到

$\sum _{d=1} ^{n} \mu(d)\left \lfloor \frac{n}{d} \right \rfloor \sum _{j=1}^{ \left \lfloor \frac{m}{d} \right \rfloor }f(j*d)$

由于$f(i)$是完全积性函数,所以$f(j*d)=f(j)*f(d)$

所以可以写成$\sum _{d=1} ^{n} \mu(d)f(d)\left \lfloor \frac{n}{d} \right \rfloor \sum _{j=1}^{ \left \lfloor \frac{m}{d} \right \rfloor }f(j)$

然后我成功的yy出了怎么线性筛$f(i)$的值

......然后不会处理$\sum _{j=1}^{ \left \lfloor \frac{m}{d} \right \rfloor }f(j)$这一部分了……我以为可以杜教筛这部分，结果发现不行

最后在空间允许的范围内打了时间复杂度为$O(n)$的24分算法

考完试和$wq$交流之后我发现......

$f(i)\Leftrightarrow[(i,k)==1]$

……当时整个人傻了

后来的后来发现自己式子的某个主流题解是一样的……

更加心里苦

我们还是来证明一下上面那个式子为啥成立吧……

设循环节长度为$L$

那么根据题给的循环节生成方式，应该有$x \% y=x*k^{L} \% y$

由于题面要求被计算的$x$和$y$互质，所以我们可以消去$x$

从而有$k^{L} \equiv 1(mod y)$

这就意味着$k$在模$y$的意义下有逆元

如果一个数在模意义下有逆元……我们回想一下$exgcd$的过程

那么就会有$(k,y)==1$

这就好说了……真是心累

然后还有一点我考试没有想到的，就是快速计算$F(n)=\sum _{i=1} ^{n} [(i,k)==1]$

由于我们知道$(a,b)==(a-k*b,b)$

所以$\sum _{i=1} ^{k} [(i,k)==1]==\sum _{i=k+1} ^{2*k} [(i,k)==1]$

由此我们可以归纳出

$F(n)=\left \lfloor \frac{n}{k} \right \rfloor F(k)+F(n \% k)$

我们回到刚才我演的那个式子

$\sum _{d=1} ^{n} \mu(d) [(k,d)==1] \left \lfloor \frac{n}{d} \right \rfloor F( \left \lfloor \frac{m}{d} \right \rfloor ) $

这样，我们预处理1～k内的$F$值，就可以$O(n)$回答询问了，可以得到84分

然而事实证明剩下的16分分值和难度不成正比……

我们考虑上面这个东西，$\left \lfloor \frac{n}{d} \right \rfloor$和$\left \lfloor \frac{m}{d} \right \rfloor$是可以除法分块的

那么我们如果能处理出来一个$\sum _{i=1} ^{n} \mu(i) f(i)$前缀和就能$O(\sqrt{n})$回答了……

这里面的f(i)是我上面定义的$f(i)\Leftrightarrow[(i,k)==1]$

我们考虑怎么处理这个前缀和……

设其为$G(n,k)$,考虑化简为子问题递归解决

对于$k$的某一个质因子$p$，我们可以把$k$写成$k=p^{a}q((p,q)==1)$的形式

那么如果$(i,k)==1$,则$(i,q)==1$且$(i,p)==1$

如果我们处理出$G(n,q)$，再减去与$q$互质但是与$p$不互质的数的个数，就可以得到$G(n,k)$了

我们假设某个与$q$互质但是与$p$不互质的数$i$可以表示为$i=p^{c}y((y,q)==1,(y,p)==1)$

由于与$p$不互质，所以$c>0$又因为$c>=2$时$\mu(i)==0$所以$c==1$

因此我们考虑的是一系列形如$i=py$的$i$

我们写一下这个式子……

$G(n,k)=G(n,q)-\sum _{y=1} ^{\left \lfloor   \frac{n}{p}   \right \rfloor } \mu(py) [(y,q)==1] [(y,p)==1]$

因为$(p,y)==1,\mu(py)=\mu(p)\mu(y)$

因此就会有$G(n,k)=G(n,q)-\sum _{y=1} ^{\left \lfloor   \frac{n}{p}   \right \rfloor } \mu(p)\mu(y) [(y,p)==1][(y,q)==1]$

我们把$\mu(p)=-1$带入，有$G(n,k)=G(n,q)+\sum _{y=1} ^{\left \lfloor   \frac{n}{p}   \right \rfloor } \mu(y) [(y,p)==1][(y,q)==1]$

又由于我们一开始利用的条件，“如果$(i,k)==1$,则$(i,q)==1$且$(i,p)==1$”

所以$[(y,p)==1][(y,q)==1]\Leftrightarrow[(y,k)==1]$

所以原式可以写成$G(n,k)=G(n,q)+\sum _{y=1} ^{\left \lfloor   \frac{n}{p}   \right \rfloor } \mu(y) [(y,k)==1]$

即$G(n,k)=G(n,q)+G(\left \lfloor  \frac{n}{p}  \right \rfloor ,k)$

这样我们就可以记忆化搜索我们的$G$值，边界条件是$n==0$时$G(0,k)=0$,$n==1$时$G(1,k)=1$；

$k==1$时$G(n,1)=\sum _{i=1} ^{n} \mu(i)$，杜教筛之。

这个部分和上面的bzoj3309一样，都是对递归的巧妙应用

这种定义一个函数求递推式来帮助思考的方法很有效！

这样这题就被解决了！的确是优美的题目啊！

不知道为什么跑的贼慢的代码：

#include <map>
#include <cstdio>
#include <cstring>
using namespace std;
#define N 20000000
#define K 2010
#define LL long long
int n,m,k;
int prime[N/10+10],tot,mu[N+10],sum[K];
bool vis[N+10];
inline int min(int a,int b){return a<b?a:b;}
inline int gcd(int a,int b){return b==0?a:gcd(b,a%b);}
inline void init()
{
    register int i,j,tmp;
    for(i=1;i<=k;++i)sum[i]=sum[i-1]+(gcd(i,k)==1);
    for(mu[1]=1,i=2;i<=N;++i)
    {
        if(!vis[i])prime[++tot]=i,mu[i]=-1;
        for(j=1;j<=tot&&(tmp=i*prime[j])<=N;++j)
        {
            vis[tmp]=1;
            if(i%prime[j]==0){mu[tmp]=0;break;}
            mu[tmp]=-mu[i];
        }
    }
    for(i=1;i<=N;++i)mu[i]+=mu[i-1];
}
int bin[100],ge;
inline void divide(int kkk)
{
    register int i;
    for(i=1;prime[i]*prime[i]<=kkk;++i)
        if(kkk%prime[i]==0)
        {
            bin[++ge]=prime[i];
            while(kkk%prime[i]==0)kkk/=prime[i];
        }
    if(kkk>1)bin[++ge]=kkk;
}
inline int getF(int x)
{
    if(x<=k)return sum[x];
    return (x/k)*sum[k]+sum[x%k];
}
map<int,int>smu;
inline int getmu(int x)
{
    if(x<=N)return mu[x];
    if(smu.count(x))return smu[x];
    register int i,last,ret=1;
    for(i=2;i<=x;i=last+1)
        last=x/(x/i),ret-=(LL)getmu(x/i)*(last-i+1);
    return smu[x]=ret;
}
map<int,int>mmp[7];
inline int getsum(int x,int layer)
{
    if(x<=1)return x;
    if(layer==0)return getmu(x);
    if(mmp[layer].count(x))return mmp[layer][x];
    return mmp[layer][x]=getsum(x,layer-1)+getsum(x/bin[layer],layer);
}
int main()
{
    scanf("%d%d%d",&n,&m,&k);
    init(),divide(k);
    register int i,to=min(n,m),last;
    LL ans=0;
    // if(n>m)n^=m,m^=n,n^=m;
    for(i=1;i<=to;i=last+1)
        last=min(n/(n/i),m/(m/i)),
        ans+=(getsum(last,ge)-getsum(i-1,ge))*(LL)getF(m/i)*(n/i);
    printf("%lld\n",ans);
}

大概我最近做的题就这些……虽然还有另外一道模拟赛的题，式子还不错，先不放出来了……

杜教筛……怎么说呢，算是对狄利克雷卷积的灵活应用吧

作为优化的工具，不管是复杂度还是思想都很优秀

逐渐有一点点数学题的手感了知道上来无脑反演了，加油咯……