杜教筛

杜教筛

前置技能树:积性函数

就是对于函数(f(x))

对于任意两个互质整数(a,b),如果有(f(a)·f(b)=f(ab))

(f)为积性函数

如果对任意(a,b)成立,(f)为完全积性函数。

前置技能树:狄利克雷卷积

狄利克雷卷积是一种运算定义。

(f*g=sumlimits_{d|n}f(d)g(frac nd))

其显然满足交换律。

前言

一般我们求积性函数有优秀的(O(n))欧拉筛

但是实际运用中,我们往往需要得出积性函数的前缀和来进行运算,而且总是有一些毒瘤出题人把数据出到(1e10)之类,此时线性筛就不够用了。

为了解决这个问题,我们就需要一种新的筛法—杜教筛。

杜教筛是一种筛法,能够以(O(n^{frac23}))的时间复杂度求积性函数的前缀和

此外好像还有(min\_25)筛,复杂度为(O(n^{frac{3/4}{log_n}}))

但是学不动了(Orz)

具体推导

这个东西都是套路。。。

如果不想看公式了其实翻到下面加粗加大地方背个板子也挺不错

(f(n))为你要筛的函数,(S(n)=sumlimits_{i=1}^{n}f(i))

我们构造两个积性函数(h,g),使得(h=g*f)

(sumlimits_{i=1}^n h(i)=sumlimits_{i=1}^n(g*f)(i))

(=sumlimits_{i=1}^n sumlimits_{d|i}g(d)*f(frac nd))

考虑套路,把(d)提出来,枚举倍数。

(=sumlimits_{d=1}^n g(d)sumlimits_{i=1}^{lfloorfrac nd floor}f(i))

(=sumlimits^{n}_{d=1}gleft( d ight) Sleft( lfloor dfrac {n}{d} floor ight))

我们发现这个里面已经有我们要求的(S(n))了,于是单独拿出来。

(=g(1)S(n)+sumlimits_{d=2}^n Sleft(lfloor dfrac {n}{d} floor ight))

我们移一下项,(Large g(1)S(n)=sumlimits_{i=1}^nh(i)-sumlimits_{d=2}^nS(lfloorfrac nd floor))

这就是杜教筛的套路柿子,于是大功告成,我们就沿着这个柿子递归就可以了

...吗?

还没完...

我们发现如果直接递归,复杂度是(O(n^{frac 34}))的,非常不优秀。

于是我们就先预处理出前(n^{frac 23})的数字,这样递归的复杂度也是(O(n^{frac 23}))的了,(复杂度我也不会证)

于是我们就成功的在(O(n^{frac 23}))的复杂度内求出了(f)的前缀和。

n=1e10~11都可以做

使用方法

我们发现性能瓶颈主要在(h)函数的前缀和上,所以我们需要构造一个(g)使得(h)的前缀和能在(O(1))时间复杂度上解决。

例题(1)(f=mu)

构造(g=I)

(mu * I = epsilon)

(epsilon)的前缀和显然很菜,就是1

所以(S(n)=1-sumlimits_{d=2}^nS(lfloorfrac nd floor))

例题(2)(f=varphi)

仍然构造(g=I)

(varphi * I=id)

所以(S(n)=frac{n(n+1)}2-sumlimits_{d=2}^nS(lfloorfrac nd floor))

卡常小技巧

我们平时写杜教筛的时候,筛出的结果显然要存下来。

但是这样要手写哈希(麻烦)或者使用(STL)(unorded\_map)(\_\_gnu\_pbds::cc\_hash\_table)(慢)

有没有两全其美的办法呢?

有!

我们发现求解的总是n的约数。

于是我们开一个两倍的(dp)数组。

(xle sqrt n)

返回(dp[x])

否则返回(dp[sqrt n+n/x])就可以了。

这样就可以跑得飞快,不用太费心卡常。

洛谷P4213 【模板】杜教筛

/*
@Date    : 2019-08-02 19:52:32
@Author  : Adscn (adscn@qq.com)
@Link    : https://www.cnblogs.com/LLCSBlog
*/
#include<bits/stdc++.h>
using namespace std;
#define IL inline
#define RG register
#define gi getint()
#define gc getchar()
#define File(a) freopen(a".in","r",stdin);freopen(a".out","w",stdout)
IL int getint()
{
	RG int xi=0;
	RG char ch=gc;
	bool f=0;
	while(ch<'0'||ch>'9')ch=='-'?f=1:f,ch=gc;
	while(ch>='0'&&ch<='9')xi=(xi<<1)+(xi<<3)+ch-48,ch=gc;
	return f?-xi:xi;
}
template<typename T>
IL void pi(T k,char ch=0)
{
	if(k<0)k=-k,putchar('-');
	if(k>=10)pi(k/10,0);
	putchar(k%10+'0');
	if(ch)putchar(ch);
}
const int MAXN=5000000;
const int N=INT_MAX;
const int SQRN=sqrt(N)+7;
typedef long long ll;
int n,sqrn;
bool npr[MAXN+7];
int pr[MAXN+7],cnt;
ll mu[MAXN+7];
unsigned long long phi[MAXN+7];
ll dpp[SQRN<<1],dpm[SQRN<<1];
//map<int,ll>dpp,dpm;
void init()
{
	npr[1]=1;
	mu[1]=phi[1]=1;
	for(int i=2;i<=MAXN;++i)
	{
		if(!npr[i])pr[++cnt]=i,mu[i]=-1,phi[i]=i-1;
		for(int j=1;j<=cnt&&1ll*pr[j]*i<=MAXN;++j)
		{
			npr[i*pr[j]]=1;
			if(i%pr[j]==0){mu[i*pr[j]]=0,phi[i*pr[j]]=phi[i]*pr[j];break;}
			mu[i*pr[j]]=-mu[i],phi[i*pr[j]]=phi[i]*(pr[j]-1);
		}
	}
	for(int i=1;i<=MAXN;++i)mu[i]+=mu[i-1],phi[i]+=phi[i-1];
}
inline ll getphi(int x)
{
	if(x<=MAXN)return phi[x];
	if(x<=SQRN&&dpp[x])return dpp[x];
	else if(x>SQRN&&dpp[N/x+SQRN])return dpp[N/x+SQRN];
	unsigned long long ans=1ull*x*(x+1ull)/2ull;
	for(int l=2,r;r<INT_MAX&&l<=x;l=r+1)r=x/(x/l),ans-=1ull*(r-l+1)*getphi(x/l);
	if(x<=SQRN)dpp[x]=ans;
	else dpp[N/x+SQRN]=ans;
	return ans;
}
inline ll getmu(int x)
{
	if(x<=MAXN)return mu[x];
	if(x<=SQRN&&dpm[x])return dpm[x];
	else if(x>SQRN&&dpm[N/x+SQRN])return dpm[N/x+SQRN];
	ll ans=1;
	for(int l=2,r;r<INT_MAX&&l<=x;l=r+1)r=x/(x/l),ans-=1ll*(r-l+1)*getmu(x/l);
	if(x<=SQRN)dpm[x]=ans;
	else dpm[N/x+SQRN]=ans;
	return ans;
}
int main(void)
{
	init();
	int T=gi;
	while(T--)
	{
		n=gi;
		memset(dpp,0,sizeof dpp);
		memset(dpm,0,sizeof dpm);
		printf("%lld %lld
",getphi(n),getmu(n));
	}
	return 0;
}
原文地址:https://www.cnblogs.com/LLCSBlog/p/11305973.html