Nemlit 的博客

Nemlit 的博客

By a konjac

杜教筛学习笔记

posted on 2019-02-11 20:13:44 | under 题解 |

原文地址

杜教筛的作用是可以快速求出积性函数的前缀和,相比于传统的线性筛,杜教筛可以在 $O(n^{2/3})$ 的时间求出积性函数的前缀和。

PS:以下所有除法运算为向下取整(不知道向下取整的LaTex)

前置芝士1:积性函数

什么是积性函数呢?若a⊥b(a,b互质),那么有 $f(ab)=f(a)f(b)$ 则称 $f(x)$ 的积性函数。

前置芝士2:狄利克雷卷积

定义两个积性函数的狄利克雷卷积为 $(f*g)(i)=\sum_{d|n}f(d)$ * $g(\frac{d}{n}) $

给出几个常见的狄利克雷卷积的结论

( $I(x)=1,id(x)=x,\epsilon(x)=[x==1]$ )

  1. $\mu*I=\epsilon$

证明: $\mu*I=\sum_{d|n}\mu(d)*I(\frac{n}{d})=\sum_{d|n}\mu(d)=[n==1]=\epsilon$

倒数第二步是莫比乌斯函数的一个结论

  1. $φ*I=id(n)$

证明: $φ*I=\sum_{d|n}φ(d)*I(\frac{n}{d})=\sum_{d|n}φ(d)=n=id(n)$

对于倒数第二步,考虑统计与 $n$ 的 $gcd$ 分别为 $1,2,...,n$ 的过程:

对于 $1-n$ 所有数,与 $n$ 的 $gcd$ 为 $1$ 的有 $φ(n)$ 个;若 $2|n$ ,则与 $n$ 的 $gcd$ 为2有 $φ(\frac{n}{2})$ 个……以此类推,总共有 $n$ 个数,故 $\sum_{d|n}φ(d)=n$

  1. $id*\mu=φ(n)$

证明: $φ(n)*I*\mu=id*\mu$ ----------------结论2可知

$φ(n)*\epsilon=id*\mu$ ------------------------------结论1可知

$φ(n)=id*\mu$ ----------------------------------把卷积拆开即可

前置芝士3:莫比乌斯反演

若 $f=g*I$ ,则有 $g=f*\mu$

我们把卷积拆开

原式为:若 $f(n)=\sum_{d|n}g(d)$ ,则有 $g(n)=\sum_{d|n}f(\frac{n}{d})*\mu(d)$

证明: $f*\mu=g*I*\mu=g*\epsilon=g$

所以 $1$ 的逆是 $\mu$

正文:杜教筛

我们构造积性函数f,g,令 $S(n)=\sum_{i=1}^{n}f(i)$

$\sum_{i=1}^{n}(f*g)(i)=\sum_{i=1}^{n}\sum_{d|n}f(d)*g(\frac{n}{d})=\sum_{i=1}^{n}g(i)*\sum_{j=1}^{[\frac{n}{i}]}f(j)=\sum_{i=1}^{n}g(i)*S(\frac{n}{i})$

移向得: $g(1)*S(n)=\sum_{i=1}^{n}(f*g)(i)-\sum_{i=2}^{n}g(i)*S(\frac{n}{i})$

$S(n)=\frac{\sum_{i=1}^{n}(f*g)(i)-\sum_{i=2}^{n}g(i)*S(\frac{n}{i})}{g(1)}$ //以后称此为杜教筛公式

我们发现,S(n)是我们所求,我们把S(n)拆成了递归形式,假设我们可以求出 $\sum_{i=1}^{n}(f*g)(i)$ , $g(x)$ ,我们就可以对 $\sum_{i=2}^{n}g(i)*S(\frac{n}{i})$ 进行整除分块,从而递归求出 $S(N)$ !

看懂了上述证明,那本题就很简单了吧。

一: $\sum_{i=1}^{n}\mu(i)$

因为 $\mu*I=\epsilon$ ,所以我们定义 $g=I$ ,带入杜教筛公式得 $S(n)=1-\sum_{i=2}^{n}S(\frac{n}{i})$ ,然后就是整除分块的模板了

二: $\sum_{i=1}^{n}φ(i)$

因为 $φ*I=id(n)$ (其实也可以用 $φ(n)=id*\mu$ ),所以我们还是定义 $g=I$ ,带入杜教筛公式得 $S(n)=\frac{(1+x)*x}{2}-\sum_{i=2}^{n}S(\frac{n}{i})$ ,然后发现和上面区别并不大

这样直接求显然是不能做到非线性的,所以我们需要先预处理出一部分的 $φ(i),\mu(i)$ ,然后对于我们处理过的答案,我们可以用一个map来存储(相当于记忆化搜索)以优化时间复杂度。

经过理论证明,当预处理的值= $n^{\frac{2}{3}}$ 时,时间复杂度为 $O(n^{\frac{2}{3}})$

注意几点:

1.本题卡常,所以不要全部开longlong,也不要用map,而用c++11的STL:unordered_map(少了map的排序,复杂度不带log,但常数仍然很大,而且本地要开c++11才能过编译)

2.注意杜教筛的公式, $\sum_{i=2}^{n}g(i)*S(\frac{n}{i})$ 是从2开始的,不然会进入死递归(我也不知道叫什么)

代码如下:

#include<bits/stdc++.h>
using namespace std;
#define il inline
#define re register
#define debug printf("Now is Line : %d\n",__LINE__)
#define file(a) freopen(#a".in","r",stdin);freopen(#a".out","w",stdout)
#define ll long long
#define mod 1000000007
il int read()
{
    re int x=0,f=1;re char c=getchar();
    while(c<'0'||c>'9') {if(c=='-') f=-1;c=getchar();}
    while(c>='0'&&c<='9') x=x*10+c-48,c=getchar();
    return x*f;
}
#define maxn 5000000
ll pai[maxn+5];
int prim[maxn+5],cnt,mu[maxn+5];
bool vis[maxn+5];
unordered_map<int,int> ans_mu;
unordered_map<int,ll> ans_pai;
il void init()
{
    mu[1]=pai[1]=1;
    for(re int i=2;i<=maxn;++i)
    {
        if(!vis[i]) prim[++cnt]=i,mu[i]=-1,pai[i]=i-1;
        for(re int j=1;j<=cnt&&prim[j]*i<=maxn;++j)
        {
            vis[prim[j]*i]=1;
            if(i%prim[j]==0)
            {
                pai[i*prim[j]]=pai[i]*prim[j];
                break;
            }
            pai[i*prim[j]]=pai[prim[j]]*pai[i];
            mu[i*prim[j]]=-mu[i];
        }
    }
    for(re int i=1;i<=maxn;++i) mu[i]+=mu[i-1],pai[i]+=pai[i-1];
}
il ll get_pai(ll x)
{
    if(x<=maxn) return pai[x];
    if(ans_pai[x]) return ans_pai[x];
    ll ans=((1ll+x)*x)/2ll;
    for(re int l=2,r;l<=x;l=r+1)//其实这里可能会爆int,可以改用用unsigned int
    {
        r=x/(x/l);
        ans-=1ll*(r-l+1)*get_pai(x/l);
    }
    return ans_pai[x]=ans;
}
il int get_mu(int x)
{
    if(x<=maxn) return mu[x];
    if(ans_mu[x]) return ans_mu[x];
    int ans=1;
    for(re int l=2,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        ans-=(r-l+1)*get_mu(x/l);
    }
    return ans_mu[x]=ans;
}
il void doit()
{
    int T=read();
    while(T--)
    {
        int x=read();
        printf("%lld %d\n",get_pai(x),get_mu(x));
    }
}
signed main()
{
    init(),doit();
    return 0;
}