[模板] 杜教筛 && bzoj3944-Sum

杜教筛

浅谈一类积性函数的前缀和 - skywalkert's space - CSDN博客

杜教筛可以在(O(n^{frac 23}))的时间复杂度内利用卷积求出一些积性函数的前缀和.

算法

给定(f(n)), 现要求(S(n)=sum_{i=1}^n f(i)).

定义卷积运算 ((f*g)(n) = sum_{d | n} f(d) g(frac{n}{d})).

如果存在(g(n)), 满足(f*g=h), 且(g)和(h)都能 (O(1)) 求出前缀和, 我们可以较快地求出(S(n)).

注意到

[egin{aligned} sumlimits_{i=1}^{n}(f*g)(i) &= sumlimits_{i=1}^{n} sum limits _{d|i} f(d)g(frac{i}{d}) \ &= sum limits _{d=1}^{n} g(d)sumlimits _{i=1}^{lfloor frac{n}{d} floor } f(i) \ &= sum limits _{d=1}^{n} g(d) S(lfloor frac{n}{d} floor) end{aligned} ]

因此, 有

[g(1)S(n)=sumlimits_{i=1}^{n}(f*g)(i) - sum limits _{i=2}^{n} g(i) S(lfloor frac{n}{i} floor) ]

可以递归(并记忆化)下去.

对于复杂度: 展开一层递归, 通过积分可以求出时间复杂度为 (O(n^{frac 34})).

如果预处理前 (m) 个答案, 利用同样的方法可以得到复杂度为 (O(m + frac n{sqrt m})), 当 (m = n^{frac 23}) 时取最小值为 (O(n^{frac 23})).

并不知道为什么算复杂度时可以只展开一层

Upd: 关于为什么算复杂度时只展开一层:
递归的 (x) 显然为 (lfloor frac ni floor) 的形式, 可以通过哈希表(或者下面的另一种方法)存储. 那么递归到第二层的时候会发现要求的值已经求过了, 因此只需展开一层就行了.

关于卷积

显然, 卷积运算满足交换律和结合律, 可以推式子验证一下.

另外, 积性函数的卷积仍然为积性函数.

定义函数 (epsilon(n) = [n=1], I(n) = 1, id(n) = n), 有

[f * epsilon = f ]

这是(epsilon)函数的定义.

[phi * I = id ]

[mu * I = epsilon ]

这是莫比乌斯函数的定义.

[mu * id = phi ]

[id * id = id cdot d ]

[(phi cdot id) * id = id^2 ]

Problem Description

求 (phi (n)) 和 (mu (n)) 的前缀和. (1 le n le 2^{31}-1).

Code

另外, 卡常数...
用long long会TLE, 改成unsigned int就不会.
似乎不少毒瘤数论题都卡常

#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<set>
#include<map>
#include<unordered_map>
using namespace std;
#define rep(i,l,r) for(register int i=(l);i<=(r);++i)
#define repdo(i,l,r) for(register int i=(l);i>=(r);--i)
#define il inline
typedef double db;
typedef long long ll;
typedef unsigned long long ull;
typedef unsigned int uint;

//---------------------------------------
const int blsz=5e6+50,sqsz=5e4+50;
ll bnd=5e6;
ll t,n;

int nopr[blsz],pr[blsz],pp=0;
ll mu[blsz],phi[blsz];
void init(){
	nopr[1]=mu[1]=phi[1]=1;//a
	rep(i,2,bnd){
		if(nopr[i]==0)pr[++pp]=i,mu[i]=-1,phi[i]=i-1;//b
		rep(j,1,pp){
			if(i*pr[j]>bnd)break;
			nopr[i*pr[j]]=1;
			if(i%pr[j])mu[i*pr[j]]=-mu[i],phi[i*pr[j]]=phi[i]*phi[pr[j]];
			else{mu[i*pr[j]]=0,phi[i*pr[j]]=phi[i]*pr[j];break;}
		}
	}
	rep(i,1,bnd)phi[i]+=phi[i-1],mu[i]+=mu[i-1];
}

//ll sq,vmu[n12sz*2],vphi[n12sz*2];
//ll tr(ll v){return v<}

unordered_map<uint,ll> ansmu,ansphi;

ll getphi(uint n){
	if(n<=bnd)return phi[n];
	ll &tmp=ansphi[n];
	if(tmp)return tmp;
	ull res=(ull)n*(n+1)/2;
	for(uint l=2,r;l<=n;l=r+1){//using unsigned int instead of ll
		r=n/(n/l);
		res-=(ll)(r-l+1)*getphi(n/l);
	}
	return tmp=res;
}

ll getmu(uint n){
	if(n<=bnd)return mu[n];
	ll &tmp=ansmu[n];
	if(tmp)return tmp;
	tmp=1;
	for(uint l=2,r;l<=n;l=r+1){
		r=n/(n/l);
		tmp-=(r-l+1)*getmu(n/l);
	}
	return tmp;
}

int main(){
//	freopen("a.in","r",stdin);
	ios::sync_with_stdio(0),cin.tie(0);
	init();
	cin>>t;
	rep(cs,1,t){
		cin>>n;
		cout<<getphi(n)<<' '<<getmu(n)<<'
';
	}
	return 0;
}

unordered_map的地方也可以换为

struct tmap{
	ll a[sqsz],b[sqsz];
	void cl(){memset(b,0,sizeof(b));}
	ll& operator[](int p){
		if(p<5e4)return a[p];
		else return b[n/p];
	}
}ansmu,ansphi;

但是并没有变快... 可能是unordered_map常数小吧...