Min_25 筛 | yijan's Blog

Min_25 筛也是一种求积性函数前缀和的算法。

Min_25 筛的应用前提是可以找到一个完全积性函数在质数处和要求的函数相同。所以我们假装得到了一个完全积性函数 $f'$ 它在质数处的取值和 $f$ 相同。

首先，我们把 $1$ 扣掉，最后把 $1$ 加回来即可。所以后面的范围的下界都是 $2$ 。

我们考虑把所求拆一下，然后在非质数部分枚举一下最小的质因子以及这个质因子的次数：

\sum_{i \in P} f(i) + \sum_{p^e \le n , p \le \sqrt n} f(p^e) \left( \sum_{i \le \lfloor\frac{n}{p^e}\rfloor,\text{minp} > p} f(i) \right)

考虑把质数部分和合数部分分开计算。进行一个 $dp$ ：

g(n,k) = \sum_{i \le n , \text{minp} > p_k \or i \in P} f'(i)

即 $2 \sim n$ 中质数或最小质因子大于 $p_k$ 的合数的 $f'$ 的和。

考虑如何用 $g(n,k-1)$ 去推 $g(n,k)$ 。不难发现我们减掉了一些数，我们去掉的是最小质因子恰好是 $p_k$ 的合数。我们考虑从这样的数里面提出来一个 $p_k$ ，然后剩下的数的 $f'$ 的和一定是

g(\lfloor \frac{n}{p_k} \rfloor,k-1) - g(p_{k-1},k-1)

这个式子中 $g(p_{k-1},k-1)$ 显然就是所有 $2 \sim p_{k-1}$ 的质数的和。因为我们提出了 $p_k$ 剩下的最小质因子肯定不能是一个 $< p_k$ 的质数，所以要把这些数减掉。

也就是说，转移就是

g(n,k) = g(n,k - 1) - f'(p_k) \left( g(\lfloor \frac{n}{p_k} \rfloor,k-1) - g(p_{k-1},k-1) \right)

正是因为 $f'$ 是一个完全积性函数，在这里可以直接提出去一个 $p_k$ ，剩下的里面可能还有 $p_k$ 。

考虑如果 $k \ge \sqrt n$ ，显然 $g(n,k)$ 就是 $2 \sim n$ 的所有质数的值的和。

但是显然我们不可能对于所有 $n$ 计算这个东西。可以发现，由于

\lfloor\frac{\lfloor\frac{n}{a}\rfloor}{b}\rfloor = \lfloor\frac n {ab}\rfloor

我们可以直接对所有 $\lfloor \frac{n}{a} \rfloor$ 的形式计算 $g(x,k)$ 。

然后考虑如何计算答案。仍然可以考虑 $dp$ ：

h(n,x) = g(n) - sp_x + \sum_{p_k^e \le n,k > x} f(p_k^e) \left( h(\lfloor\frac{n}{p_k^e}\rfloor,k) + [e \neq 1] \right)

由于质数被提前了，所以 $k$ 只用枚举到 $O(\sqrt n)$ 。

总复杂度是 $O(\frac{n^{\frac 3 4}}{\log n})$ ，我也不知道怎么分析出来的。

一些实现细节

首先，在质数部分需要保存 $\lfloor \frac n x\rfloor$ 。

如果 $x \le \sqrt n$ ，我们可以直接把 $\lfloor \frac n x \rfloor$ 的值放到 $A1[x]$ 里。

如果 $x > \sqrt n$ ，我们可以把 $\lfloor \frac n x\rfloor$ 放到 $A2[\lfloor \frac n x\rfloor]$ 里。

然后 $g$ 显然可以滚动数组，因为真正有用的只有 $g(n,\sqrt n)$ 这个东西。

然后由于某些原因 Min_25 筛中是没有必要进行记忆化的。

#include "iostream"
#include "algorithm"
#include "cstring"
#include "cstdio"
#include "cmath"
#include "vector"
#include "map"
#include "set"
#include "queue"
using namespace std;
#define MAXN 200006
//#define int long long
#define rep(i, a, b) for (int i = (a), i##end = (b); i <= i##end; ++i)
#define per(i, a, b) for (int i = (a), i##end = (b); i >= i##end; --i)
#define pii pair<int,int>
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define eb emplace_back
#define vi vector<int>
#define all(x) (x).begin() , (x).end()
#define mem( a ) memset( a , 0 , sizeof a )
typedef long long ll;
const int P = 1e9 + 7;
const int iv2 = P + 1 >> 1 , iv3 = ( P + 1 ) / 3;
ll n;

namespace Min_25 {
	ll n; int B;
	int pri[MAXN] , en , sp[MAXN] , sp2[MAXN];
	void sieve( ) {
		B = sqrt( n + 0.5 );
		rep( i , 2 , B ) {
			if( !pri[i] ) pri[++ en] = i , sp[en] = ( sp[en - 1] + i ) % P , sp2[en] = ( sp2[en - 1] + i * 1ll * i ) % P;
			for( int j = 1 ; j <= en && i * pri[j] <= B ; ++ j ) {
				pri[i * pri[j]] = 1;
				if( i % pri[j] == 0 ) break; 
			}
		}
	}
	int cc( ll x ) { x %= P; return ( x * ( x + 1 ) / 2 + P - 1 ) % P; }
	int cc2( ll x ) { x %= P; return ( x * ( x + 1 ) / 2 % P * ( 2 * x + 1 ) % P * iv3 + P - 1 ) % P; }
	ll g1[MAXN] , g2[MAXN] , A[MAXN];
	int A1[MAXN] , A2[MAXN];
	int po( ll x ) {
		return x <= B ? A1[x] : A2[n / x];
	}
	void getG( ) {
		int tot = 0;
		for( ll l = 1 , r ; l <= n ; l = r + 1 ) {
			r = n / ( n / l );
			++ tot;
			A[tot] = n / l;
			if( A[tot] <= B ) A1[A[tot]] = tot; else A2[n / A[tot]] = tot;
			g1[tot] = cc( A[tot] ) , g2[tot] = cc2( A[tot] );
		}
		rep( i , 1 , en )
			for( int j = 1 ; pri[i] * 1ll * pri[i] <= A[j] ; ++ j ) {
				g1[j] = ( g1[j] + P - pri[i] * 1ll * ( g1[po( A[j] / pri[i] )] + P - sp[i - 1] ) % P ) % P;
				g2[j] = ( g2[j] + P - pri[i] * 1ll * pri[i] % P * ( g2[po( A[j] / pri[i] )] + P - sp2[i - 1] ) % P ) % P;
			}
	}
	
	int g( int x ) {
		return ( g2[x] < g1[x] ? g2[x] - g1[x] + P : g2[x] - g1[x] );
	}
	int f( ll x ) {
		x %= P;
		return x * ( x - 1 ) % P;
	}
	
	int S( ll n , int k ) {
		if( n <= pri[k] ) return 0;
		int re = ( 1ll * g( po( n ) ) + P - sp2[k] + sp[k] ) % P;
		for( int i = k + 1 ; pri[i] * 1ll * pri[i] <= n && i <= en ; ++ i ) {
			ll cur = pri[i];
			for( int e = 1 ; cur <= n ; ++ e ) {
				re = ( re + f( cur ) * 1ll * ( S( n / cur , i ) + ( e != 1 ) ) ) % P;
				cur *= pri[i];
			}
		}
		return re;	
	}
	
	int solve( ll x ) {
		n = x;
		sieve( );
		getG( );
		return ( S( n , 0 ) + 1 ) % P;
	}
}

void solve() {
	cin >> n;
	cout << Min_25::solve( n ) << endl;
}

signed main() {
//    int T;cin >> T;while( T-- ) solve();
    solve();
}