用杜教筛求解数论函数前缀和

杜教筛用来求数论函数(f)前缀和。复杂度为(O(n^{frac{2}{3}}))

前提

如果我们要求(S(n)=sumlimits_{i=1}^nf(i))，那么需要找到一个数论函数(g)，满足(g)的前缀和可以非常快速的求出来，并且(g*f)的前缀和可以非常快速的求出来。

推导

既然(g*f)的前缀和可以非常快速的求出来，我们就求(g*f)的前缀和。

即(sumlimits_{i=1}^nsumlimits_{d|i}g(frac{i}{d})f(d))。

然后我们想得到的是(sumlimits_{i=1}^nf(i))。所以我们让上面的式子减去一个(sumlimits_{i=1}^nsumlimits_{d|i,d eq i}g(frac{i}{d})f(d))。

对于后面这个式子,我们用(frac{n}{d})来代替(d)，就变成了

[sumlimits_{i=1}^nsumlimits_{d eq 1,d|i}g(d)f(frac{i}{d}) ]

[=sumlimits_{d=2}^n g(d)sumlimits_{i=1}^{lfloorfrac{n}{d} floor}f(i) ]

[=sumlimits_{d=2}^n g(d)S(lfloorfrac{n}{d} floor) ]

因为(g)的前缀和可以快速求出，所以直接数论分块，后面的(S(frac{n}{d}))直接递归就好了。

这样我们得到的是(sumlimits_{i=1}^ng(1)f(i)=g(1)sumlimits_{i=1}^n f(i))，所以答案除以(g(1))（一般为1）就好了。

例子

以求(varphi)的前缀和为例。因为(f*1=Id)，(1)和(Id)的前缀和都非常好求，所以我们令(g)为(1)即可。

[S(n)=sumlimits_{i=1}^nsumlimits_{d|i}varphi(d)-sumlimits_{i=1}^nsumlimits_{d|i,d eq i}varphi(d)\ = sumlimits_{i=1}^ni-sumlimits_{i=1}^nsumlimits_{d|i,d eq 1}varphi(frac{n}{d})\ =frac{n(n+1)}{2}-sumlimits_{d=2}^nsumlimits_{i=1}^{lfloorfrac{n}{d} floor}varphi(i)\ = frac{n(n+1}{2}-sumlimits_{d=2}^nS(lfloorfrac{n}{d} floor) ]

再来推一下(mu)的前缀和。因为(mu * 1= epsilon)，(1)和(epsilon)的前缀和都非常好求，所以还是令(g=1)

[S(n)=sumlimits_{i=1}^nsumlimits_{d|i}mu(d)-sumlimits_{i=1}^nsumlimits_{d|i,d eq i}mu(d)\ =1-sumlimits_{d=2}^nsumlimits_{i=1}^{lfloorfrac{n}{d} floor}mu(i)\ = 1-sumlimits_{d=2}^nS(lfloorfrac{n}{d} floor) ]

关于预处理

这样直接搜的复杂度是(O(n^{frac{3}{4}}))，为了使复杂度更优，我们需要先预处理出一部分答案，如果我们预处理除了([1,K])的答案，当计算([1,K])中的结果时，直接返回即可。

可以证明当(K)取(n^{frac{2}{3}})时，复杂度最优秀为(n^{frac{2}{3}})

关于记忆化

为了让复杂度是正确的，我们肯定要将每次算出的结果记忆化下来。因为(n)比较大，所以需要用(map)来记忆化。这样复杂度就会多个(log)

还有一种方法，因为我们每次递归到的(x)肯定满足:所有满足(frac{n}{y}=frac{n}{x})的(y)中，只有(x)会被计算到。所以我们可以用一个数组ma记忆化，当查询的(n)小于等于(K)时，我们直接范围答案，当查询的(n)大于(K)时，我们查看(ma[lfloor frac{n}{K} floor])中的值即可。

当有多次询问时，第二种方法需要清空，复杂度可能不如第一种。

代码

以求(varphi)的前缀和为例。

void pre() {
	phi[1] = 1;
	for(int i = 2;i < N;++i) {
		if(!vis[i]) {
			prime[++tot] = i;
			phi[i] = i - 1;
		}
		for(int j = 1;j <= tot && prime[j] * i < N;++j) {
			vis[prime[j] * i] = 1;
			if(i % prime[j]) {
				phi[i * prime[j]] = phi[i] * (prime[j] - 1);
			}
			else {
				phi[i * prime[j]] = phi[i] * prime[j];
				break;
			}
		}
	}
	for(int i = 2;i < N;++i) phi[i] = (phi[i] + phi[i - 1]) % mod;
}
ll MAX;

ll PHI(ll n) {
	if(n < N) return phi[n];
	
	if(vis[MAX / n]) return maphi[MAX / n];
	vis[MAX / n] = 1;
	ll ret = (n % mod) * ((n + 1) % mod) % mod * inv % mod;

	for(ll l = 2,r;l <= n;l = r + 1) {
		r = n / (n / l);
		ret -= ((r - l + 1) % mod) * PHI(n / l) % mod;
		ret %= mod;
	}
	return maphi[MAX / n] = ret;
}