「学习笔记」Min25筛
前言
周指导今天模拟赛五分钟秒第一题,十分钟说第二题是 ( ext{Min25}) 筛板子题,要不是第三题出题人数据范围给错了,周指导十五分钟就 ( ext{AK}) 了,为了向 ( ext{AK})王 学习,真诚的膜拜他,接受红太阳的指导,下午就学习了一下 ( ext{Min25}) 筛。
简介
如果 (f(n)) 是一个积性函数,且 (f(n)) 是一个关于 (n) 的简单多项式,并可以快速算出 (f(p^k), p ext{ is prime, } k geq 0) 的值,那么大概可以用 ( ext{Min25}) 筛在 (mathcal O(frac{n^frac{3}{4}}{log_n})) 求出 (sum_{i=1}^nf(i)) 的值。
算法
首先令 (s) 为所有 (iin[1,n], lfloorfrac{n}{i} floor) 的集合,并试图预处理出 (forall x in s sum_{i=2}^x [i ext{ is prime}] f(i))
可能当 (x) 不是质数的时候,(f(x)) 不太好求,我们先假装所有数都是质数,然后用类似埃式筛法的过程把不是质数的 (f(x)) 值给筛掉。
具体的算法之前,先定义一些东西:
令 (P={prime_1dots prime_m}) 表示前 (m) 个质数的集合,并且 (prime_{m+1}>sqrt n) 。
令 $g(n,k) ={sum_{i=2}^n[i ext{ is prime} ext{ or } minp(i)>prime_k]}f(i) $ 表示前 (n) 个数中满足是质数或者最小质因子大于第 (k) 个质数的数假装它为质数求出来的 (f(x)) 之和。
(g(n, k)) 的本质是筛掉了前 (n) 个数筛去前 (k) 个质数的倍数后剩下的那些数的 (f(x)) 之和,显然,(g(S,m)) 就是要预处理的东西。
可以通过这个东西求 (g) 了
第一个转移显然,(prime_k) 筛不掉任何数了,第二个转移考虑把所有在前 (k-1) 轮最小质因子不为 (prime_k) 的数贡献减掉,因为 (f(x)) 是积性的直接搞就好了,因为 (g(n,k)) 仍然要保留前 (k-1) 个质数的贡献,所以还要加回来。
现在已经预处理出了 (forall x in s sum_{i=2}^x [i ext{ is prime}] f(i)) 的值,也就是我们求出了质数的答案,接下来通过从小到大加入质因子来求出合数的答案。
令 (S(n,k)=sum_{i=2}^n [minp(i)geq k] f(i)) ,这里的 (f(i)) 就是真实的值了,不假装他是质数,那么答案就是 (f(1)+S(n,1)) 。
质数的贡献之前已经算出来了,就是 (g(n,m)-sum_{i=1}^{k-1}f(prime_i))。
考虑枚举每个数的最小质因子以及它们幂次得到合数的转移
其中 (S(lfloorfrac{n}{prime_j^t} floor,j+1) imes f(prime_{j}^t)) 算的是后面还有别的质因子的情况,(f(prime_j^{t+1})) 算的是后面没有别的质因子的情况,可以感性理解一下。
例题
「Loj6235」 区间素数个数
求 ([1,n]) 的素数个数,(n leq10^{11})
令 (f(x)) 当 (x) 是质数的时候为 (1) ,(x) 是其它数的时候为 (0) ,虽然不满足这个东西是一个关于 (x) 的简单多项式,但是只求 (f) 是质数的时候的答案显然是对的,(g(n,m)) 就是答案。
code
/*program by mangoyang*/
#include <bits/stdc++.h>
#define inf (0x7f7f7f7f)
#define Max(a, b) ((a) > (b) ? (a) : (b))
#define Min(a, b) ((a) < (b) ? (a) : (b))
typedef long long ll;
using namespace std;
template <class T>
inline void read(T &x){
int ch = 0, f = 0; x = 0;
for(; !isdigit(ch); ch = getchar()) if(ch == '-') f = 1;
for(; isdigit(ch); ch = getchar()) x = x * 10 + ch - 48;
if(f) x = -x;
}
#define ID(x) ((x) <= T ? id[0][(x)] : id[1][n/(x)])
const int N = 1000005;
ll prime[N], id[2][N], a[N], g[N], n, m, tot, T;
int b[N];
int main(){
read(n), T = sqrt(n);
for(ll l = 1; l <= n; l = n / (n / l) + 1){
a[++m] = n / l, g[m] = a[m] - 1;
id[a[m]<=T?0:1][a[m]<=T?a[m]:n/a[m]] = m;
}
for(int i = 2; i <= T; i++){
if(!b[i]) prime[++tot] = i;
for(int j = 1; j <= tot && i * prime[j] <= T; j++){
b[i*prime[j]] = 1;
if(i % prime[j] == 0) break;
}
}
for(int i = 1; i <= tot; i++)
for(int j = 1; j <= m && prime[i] * prime[i] <= a[j]; j++)
g[j] -= g[ID(a[j]/prime[i])] - (i - 1);
cout << g[ID(n)] << endl;
return 0;
}
「51NOD1222」 最小公倍数计数
令 (f(n)=sum_{i=1}^nsum_{j=i}^n [ ext{lcm}(i,j)=n]),求 (sum_{i=a}^b f(i),aleq bleq 10^{11}) 。
首先先令 (f'(n)=sum_{i=1}^nsum_{j=1}^n[ ext{lcm}(i,j)=n]) ,显然 (f(n)=frac{f'(n)+n}{2}) 。
最后一步考虑组合意义,每一种质因子,要么全部归 (i) 要么全部归 (j) ,然后你会发现,(f'(n)) 是一个积性函数,且满足 (f(prime_i^k)=2k+1) ,直接 ( ext{Min25}) 筛即可。
code
#include<bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N = 1000005;
namespace Min25{
#define id(x) ((x) <= T ? id1[x] : id2[n/(x)])
ll prime[N];
int b[N], tot, id1[N], id2[N], m;
ll a[N], g[N], T, n;
inline void init(){
m = tot = 0, T = sqrt(n + 0.5);
for(int i = 2; i <= T; i++){
if(!b[i]) prime[++tot] = i;
for(int j = 1; j <= tot && i * prime[j] <= T; j++){
b[i*prime[j]] = 1;
if(i % prime[j] == 0) break;
}
}
for(ll l = 1; l <= n; l = n / (n / l) + 1){
a[++m] = n / l;
if(a[m] <= T) id1[a[m]] = m; else id2[n/a[m]] = m;
g[m] = 3ll * (a[m] - 1);
}
for(int j = 1; j <= tot; j++)
for(int i = 1; i <= m && prime[j] * prime[j] <= a[i]; i++){
g[i] -= g[id(a[i]/prime[j])] - 3ll * (j - 1);
}
}
inline ll solve(ll a, ll b){
if(a < prime[b]) return 0;
ll ans = g[id(a)] - 3 * (b - 1);
for(int j = b; j <= tot && prime[j] * prime[j] <= a; j++)
for(ll p = prime[j], t = 1; p * prime[j] <= a; t++, p = p * prime[j])
ans += solve(a / p, j + 1) * (2 * t + 1) + 2 * t + 3;
return ans;
}
inline ll gao(ll x){
if(x == 0) return 0;
if(x == 1) return 1;
n = x, init();
return (solve(n, 1) + 1 + n) / 2ll;
}
}
int main(){
ll a, b;
cin >> a >> b;
cout << Min25::gao(b) - Min25::gao(a - 1) << endl;
return 0;
}