杜教筛

标签：卷积   ==   clu   lld   break   min   using   reg   spl   

神仙的算法

我们如果要求

\[\sum_{i=1}^N\mu(i)\]

应该怎么办

线筛显然是最常规的操作了，但是复杂度是\(O(N)\)的，如果大一点就挂了

这个时候就需要杜教筛这种神奇的东西了，可以在非线性时间内求积性函数的前缀和

比如说我们要求的是\(f\)吧

我们设一个函数\(g\)，同时还有\(h=f\times g\)

我们设\(S(i)=\sum_{k=1}^if(k)\)

于是我们要求的就是\(S(N)\)

根据卷积的定义

\[\sum_{i=1}^Nh(i)=\sum_{i=1}^N\sum_{d|i}f(\frac{i}{d})g(d)\]

常规套路把\(d\)放到外面来枚举

\[\sum_{i=1}^Nh(i)=\sum_{d=1}^Ng(d)\sum_{d|i}f(\frac{i}{d})=\sum_{d=1}^Ng(d)\sum_{i=1}^{\frac{N}{d}}f(i)=\sum_{d=1}^Ng(d)S(\left \lfloor \frac{N}{d} \right \rfloor)\]

我们把\(d=1\)的情况单独拿出来

\[\sum_{i=1}^Nh(i)=S(n)g(1)+\sum_{d=2}^Ng(d)S(\left \lfloor \frac{N}{d} \right \rfloor)\]

于是就有

\[S(N)g(1)=\sum_{i=1}^Nh(i)-\sum_{d=2}^Ng(d)S(\left \lfloor \frac{N}{d} \right \rfloor)\]

但是这个样子有什么用呢，显然没什么用啊

但是如果你的\(g,h\)是一些非常好求的函数呢，比如说\(I,id,ε\)

那就可以为所欲为了

比如说对\(\mu\)求前缀和

\[\mu\times I=ε\]

我们把上面那一套东西带进去

\[S(N)=\sum_{i=1}^Nε(i)-\sum_{d=2}^NI(d)S(\left \lfloor \frac{N}{d} \right \rfloor)\]

我们发现这是一个递归定义的柿子啊，我们可以递归去求\(S(\left \lfloor \frac{N}{d} \right \rfloor)\)啊，之后用整除分块得到\(S(N)\)的值

自然\(\varphi\)也可以这样求前缀和

\[\varphi\times I=id\]

于是

\[S(N)=\sum_{i=1}^Nid(i)-\sum_{d=2}^NI(d)S(\left \lfloor \frac{N}{d} \right \rfloor)\]

一个套路

下面是筛\(\mu,\varphi\)的板子

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<tr1/unordered_map>
#define re register
#define maxn 5000005
#define LL long long
using namespace std::tr1;
unordered_map<int,int> ma;
unordered_map<int,LL> Ma;
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
inline int read()
{
    char c=getchar();
    int x=0;
    while(c<‘0‘||c>‘9‘) c=getchar();
    while(c>=‘0‘&&c<=‘9‘) x=(x<<3)+(x<<1)+c-48,c=getchar();
    return x;
}
int f[maxn],p[maxn>>1],mu[maxn];
LL phi[maxn];
int N[101];
int M,T;
int solve(int x)
{
    if(x<=M) return mu[x];
    if(ma.find(x)!=ma.end()) return ma[x];
    int now=1;
    for(re int l=2,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        now-=(r-l+1)*solve(x/l);
    }
    return ma[x]=now;
}
LL work(int x)
{
    if(x<=M) return phi[x];
    if(Ma.find(x)!=Ma.end()) return Ma[x];
    LL ans=(LL)x*(x+1)/2;
    for(re int l=2,r;l<=x;l=r+1)
    {
        r=x/(x/l);
        ans-=(LL)(r-l+1)*work(x/l);
    }
    return Ma[x]=ans;
}
int main()
{
    T=read();
    for(re int i=1;i<=T;i++) N[i]=read(),M=max(M,N[i]);
    M=5000000;
    f[1]=mu[1]=phi[1]=1;
    for(re int i=2;i<=M;i++)
    {
        if(!f[i]) p[++p[0]]=i,mu[i]=-1,phi[i]=i-1;
        for(re int j=1;j<=p[0]&&p[j]*i<=M;j++)
        {
            f[p[j]*i]=1;
            if(i%p[j]==0)
            {
                phi[p[j]*i]=p[j]*phi[i];
                break;
            }
            mu[i*p[j]]=-1*mu[i];
            phi[i*p[j]]=phi[i]*phi[p[j]];
        }
    }
    for(re int i=1;i<=M;i++) mu[i]+=mu[i-1],phi[i]+=phi[i-1];
    for(re int t=1;t<=T;t++) printf("%lld %d\n",work(N[t]),solve(N[t]));
    return 0;
}

杜教筛

标签：卷积   ==   clu   lld   break   min   using   reg   spl   

原文地址：https://www.cnblogs.com/asuldb/p/10205660.html