Min_25 筛

简介

Min_25 可以快速解决一类积性函数的前缀和问题。

概述

求积性函数 $f(x)$ 的前缀和。当然，在《从掌握到精通》中也有提到，其实可以拓展到部分非积性函数。

总的来说，算法第一步，求出质数 $p$ 对应的 $f(p)$ 的前缀和。第二步，求出 $f(x)$ 的前缀和。从某种意义上说，这是因为质数的幂的数量很小，可以暴搜合数的质因数分解；但质数处的取值就很难搞。

Step One

一般的积性函数是复杂的。所以我们需要用一个 完全积性函数 去 拟合 原函数。只需要让 $f(p)$ 保持不变即可，因为 $\text{step one}$ 只需要求出 $\sum_{p\in\text{Prime}}f(p)$ 。
用 $g(i,j)$ 表示 $[2,i]$ 去除前 $j$ 个质数的倍数之后，剩余数字的 $f(x)$ 之和。注意：要求 $g(i,j)$ 时刻包含 $[2,i]$ 之间所有质数的 $f(x)$，也就是说，去除的是质数的 “真倍数”（不包含质数本身）。我们的目标是求出 $g(n,+\infty)$，也就是只留下了质数。

设第 $x$ 大的质数是 $p_x$。转移过程就是去掉 $kp_j(k>1)$即可。显然 $k$ 不应该是前 $(j{\rm-}1)$ 个质数之一的倍数（含自身），因为那是已经被移除的数。又因为 $f(x)$ 是完全积性函数，有 $f(kp_j)=f(p_j)f(k)$,于是有转移式

$$g(i,j)=g(i,j−1)−f(p j)[g(⌊ i/p_j⌋,j−1)−sumf(j−1)](p j ⩽ √i)$$

这里 $\operatorname{sumf}(j{\rm-}1)=\sum_{i=1}^{j-1}f(p_i)$ ，因为 $g$ 里面没有挖掉质数本身，减去它以修正。

为何有范围限制呢？如果 $i<p_j^{\thinspace 2}$ ，那么显然没有更多的数字需要被筛掉，此时 $g(i,j)=g(i,j{\rm-}1)$ 。若在 $j$ 这一维上做滚动数组，则无需修改这些位置。同时可知，需要的质数最大是 $\sqrt n$，所以直接线性筛即可求出 $\operatorname{sumf}$ 了。

我们要对 $g(i,0)$ 赋初值。所以又涉及一个问题是，拟合函数的前缀和要易于求解。一般来说，我们会选择单项式（如果是多项式，就单项式分别做 $\text{step one}$ 然后相加）作为拟合函数。

而 $i$ 的范围仍然是 $[1,n]$，这可不行。注意到我们只需要$g(n,+\infty)$，所以 只会递归到形如 $i=\lfloor{n\over v}\rfloor$ 的地方。只对这样的 $i$ 求解 $g$ 就行啦！具体实现的时候，可以先预处理出所有这样的 $i$，放在数组中；无需 $\tt map$ 从值映射到下标，因为这些值的特点是 $\lfloor{n\over x}\rfloor$ 互不相同，所以 $x>\sqrt{n}$时就放在某个长度为 $\sqrt n$ 的静态数组里即可（用 $\lfloor{n\over x}\rfloor$作为下标）。

Step Two

注意：$f(1)$ 总是会被算漏。记得最后加上。

递归法

我们前面只求出了质数的求和是 $g(n,+\infty)$，所以现在我们只需要考虑合数了。

好消息：合数的最小质因数是不超过 $\sqrt{n}$的！于是咱可以枚举这玩意儿。记 $x$ 的最小质因数为 $\gamma(x)$。

记 $h(i,j)=\sum_{x=2}^{i}[\gamma(x)\geqslant p_j]\;f(x)$。那么答案就是 $h(n,1)$，求出它就完事儿了！

枚举最小质因子与其幂次，可以写出转移式

$h(i,j)=g(i,+∞)−sumf(j−1)+⌊∑​x⩾j⌋⌊p_x^k​⩽i∑k⩾1⌋​f(p_x^k​)[h(⌊i/p_x^k⌋,x+1)+[k​≠1]]$

用 $g(i,+\infty)-\operatorname{sumf}(j{\rm-}1)$ 计算了 $[2,i]$ 中的质数 $p_c\;(c\geqslant j)$，用 $f(p_x^{\thinspace k})\cdot h(\dots)$ 计算了含至少两个质因子的合数，用 $[k\ne 1]f(p_x^{\thinspace k})$ 计算了只含一个质因子的合数，可谓不重不漏。

有一个小细节：$\lfloor{i\over p_x^{\thinspace k}}\rfloor<p_x$时，$h(\dots)$ 这一项是 $0$，可以略去。所以递归时要求 $p_x^{\thinspace k}\leqslant\lfloor{i\over p_x}\rfloor$。而这个范围恰好使得 $p_x^{k+1}$就是所有只含一个质因子的合数！于是代码就比较简单。

同理，$p_x\leqslant\sqrt{i}$ 才可以转移。有趣的是，我们可以直接 递归，而且 不进行记忆化。据称，在 $\tt min25$ 可接受范围内，这种方法优于 $\Theta \left(n^{0.75}\over\ln n\right)$。详见文末。

递推法

显然我们也有

$$h(i,j)=h(i,j+1)+[∑1<p_j^t]​⩽i​f(p_j^t​)[h(⌊i/p_j^t​​⌋,j+1)+1]$$

初值应当为 $h(i,q)=g(i,+\infty)-\operatorname{sumf}(q{\rm-}1)\;(q\geqslant r)$，其中 $p_r$是大于 $\sqrt{i}$的最小质数。反正最后都要计算，干脆就保持 $p_r$为大于 $\sqrt{n}$的最小质数进行赋值吧。

处理则稍微麻烦些：对于 $p_j\leqslant\sqrt{i}$，正常转移。对于 $p_j>\sqrt{i}$，可以发现右侧的 $h$ 函数值为 $0$，提供贡献的只有那个 $f(p_j)$了，也就是较大的质数处的 $f$ 值；用 $\text{step one}$ 中求出的 $g(i,+\infty)$ 算一下就好。

这样的处理很麻烦。不如学习一下 $\text{step one}$ 的思想，令 $h(i,j)$ 始终包含全体质数！那么转移式就应该写为

$$h(i,j)=h(i,j+1)+[∑1<p_j^t]​⩽i​f(p_j^t​)[h(⌊i/p_j^t​​⌋,j+1)-sumf(j)+[t\ne1]]$$

那么 $p_j>\sqrt{i}$时已经无需转移了！初值则设置为 $h(i,+\infty)=g(i,+\infty$，相当于在 $g$ 的基础上接着递推，结束。

【模板】Min_25筛

定义积性函数$f(x)$，且$f(p^k)=p^k(p^k-1)$（$p$是一个质数），求
$\sum_{i=1}^n f(i)$
对$10^9+7$取模。
对于$100\%$的数据，保证$1\le n\le 10^{10}$

代码

时空复杂度 $O(n^{0.75}/$In$n)$

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cctype>
using namespace std;
# define rep(i,a,b) for(int i=(a); i<=(b); ++i)
# define drep(i,a,b) for(int i=(a); i>=(b); --i)
typedef long long llong;
inline int readint(){
	int a = 0, c = getchar(), f = 1;
	for(; !isdigit(c); c=getchar())
		if(c == '-') f = -f;
	for(; isdigit(c); c=getchar())
		a = (a<<3)+(a<<1)+(c^48);
	return a*f;
}

const int MOD = 1e9+7, SQRTN = 100005;
int primes[SQRTN], primes_size;
bool isPrime[SQRTN];
struct Node{
	int one, two;
	Node() = default;
	Node(int x):one(x),two(int(llong(x)*x%MOD)){}
	Node(int _o,int _t):one(_o),two(_t){}
	Node operator * (const Node &t) const {
		return Node(int(llong(one)*t.one%MOD),
			int(llong(two)*t.two%MOD));
	}
	Node operator - (const Node &t) const {
		return Node((one-t.one+MOD)%MOD,(two-t.two+MOD)%MOD);
	}
	int val() const { return (two-one+MOD)%MOD; }
};
Node sumf[SQRTN];
void sieve(int n){
	memset(isPrime+2,true,n-1);
	for(int i=2; i<=n; ++i){
		if(isPrime[i]){
			primes[++ primes_size] = i;
			Node &v = sumf[primes_size] = sumf[primes_size-1];
			if((v.one += i) >= MOD) v.one -= MOD;
			v.two = int((v.two+llong(i)*i)%MOD);
		}
		for(int j=1; j<=primes_size&&primes[j]<=n/i; ++j){
			isPrime[i*primes[j]] = false;
			if(i%primes[j] == 0) break;
		}
	}
}

const int inv2 = (MOD+1)>>1, inv3 = (MOD+1)/3;
int haxi[2][SQRTN]; // value or divisor
inline int& index_(const llong &x,const llong &n){
	return (x < SQRTN) ? haxi[0][x] : haxi[1][n/x];
}
llong w[SQRTN<<1]; ///< positions to get sum
Node g[SQRTN<<1]; ///< 2 times SQRTN
void step_one(const llong &n){
	int tot = 0; ///< allocate index
	for(llong i=1; i<=n; i=n/(n/i)+1){
		w[++ tot] = n/i; index_(n/i,n) = tot;
		g[tot].one = int(((w[tot]%MOD+1)*(w[tot]%MOD)>>1)%MOD)-1;
		g[tot].two = int((w[tot]<<1|1)%MOD*(w[tot]%MOD+1)
			%MOD*(w[tot]%MOD)%MOD*inv2%MOD*inv3%MOD)-1;
	}
	for(int j=1; j<=primes_size; ++j){
		if(primes[j] > n/primes[j]) break;
		for(int i=1; i<=tot; ++i){
			if(primes[j] > w[i]/primes[j]) break;
			g[i] = g[i]-Node(primes[j])*(
				g[index_(w[i]/primes[j],n)]-sumf[j-1]);
		}
	}
}

inline llong func(const llong &x){
	return (x-1)%MOD*(x%MOD)%MOD; // definition
}
int step_two(const llong &x,int i,const llong &n){
	if(primes[i] > x) return 0;
	int res = g[index_(x,n)].val()-sumf[i-1].val();
	for(; i<=primes_size&&primes[i]<=x/primes[i]; ++i){
		llong t = primes[i], fk = x/t;
		for(; t<=fk; t*=primes[i])
			res = int((func(t*primes[i])+func(t)
				*step_two(x/t,i+1,n)+res)%MOD);
	}
	return (res >= 0) ? res : (res+MOD);
}

int main(){
	sieve(SQRTN-1); // just do it
	llong n; scanf("%lld",&n); step_one(n);
	printf("%d\n",step_two(n,1,n)+1);
	return 0;
}

更多内容

关于 Min_25 已经有更快的算法了，详见此博客。