更新日志
2023-05-31
博客第一版
2023-06-23
更新了博客中个别地方的表达。
对“例题导入”部分作了补充。
(想加几道例题,但快生地会考了,不得不先把注意力放在生地上 qwq,只能先咕咕咕了)
2023-07-01
生地考完了,期末考完了,放一天假(今天)然后就回学校竞赛。
根据大佬 RP_INT_MAX 的要求,暂时先加一道例题,随便简单补充了下
(等到以后轻松一点再考虑要不要再加例题)
添加了目录
目录
0x00 前言
0x01 符号定义
0x02 min25 筛的用处
0x10 第一层:求解
0x11 第二层:求解
0x12 回到第一层:求解
0x20 例题导入
0x21 补充例题
0x30 结论汇总
0x40 参考资料
0x00 前言
由于本蒟蒻的搜索能力较差,很多地方在学习是没有找到明确的推导过程或定义的讲解,所以只好手推或自行定义,如有错误,欢迎指正。
本篇主要是在讲 min25 筛的各条公式如何无脑的用推式子的方式推出,所以会有很多推式子的内容,我会在结论前标注以供公式恐惧症患者快速查看。
0x01 符号定义
下文中意义不明的符号,都可以在这找到。
符号/函数 | 含义 |
---|---|
质数集合 | |
第 |
|
质数集合的子集,包含 |
|
0x02 min25 筛的用处
min25 筛被用于处理一类积性函数的前缀和。
也即,对于一个函数
一、函数
二、若
三、若
0x10 第一层:求解
如果你认真地阅读了 min25 筛的用处,你会发现使用条件几乎都是在限制质数,这是因为在 min25 筛的处理中需要用到大量的质数的性质。
我们可以将上式按照代入的自变量划分成三个部分:质数,合数与 1。
合数想与质数挂上钩,最好的方法就是分解质因数。
我们可以枚举一个合数的最小质因子及其次数。
由于
第二个求和符号枚举
上式中,为了不出现
我们不妨设
根据上面的式子,我们可以依照相同的思路得到
一、分质数与合数(这里定义的
二、合数分解质因数
还可以再展开一下(结论一
):
至此,我们就得到了答案计算的递归式了。
接下来就需要处理前面的质数部分了。
0x11 第二层:求解
显然,1 不是质数。
所以原式可以变成:
为了解决它,我们引入一个函数:
那么,
由于
可以举个例子(究竟如何线性组合需要依照具体题目进行分析):
它由
于是,我们对
接下来考虑递推
显然根据定义有:
如果
如果
否则继续讨论。
更改求和指标,使枚举
只要
将后面向
显然上式是有错的,对于大于
我们为了方便,再加入一个定义:
他经过线性组合,得到
于是:
总结一下(结论二
):
通过线性组合
0x12 回到第一层:求解
我们此前已经推出的结论:
带入上面求出的结果(最终结论
):
自此,我们已经完成了大部分的理论推导,接下来我会结合例题进行讲解。
0x20 例题导入
例题解析
请先读题:传送门
首先,他肯定适用 min25 筛。
我们可以简单的改动一下题中的式子(展开为多项式):
显然他可以由
于是,我们可以分别求出
分析完题目,我们开始做:
第一步是要处理出
根据定义:
于是:
能这样拆是因为
容易发现,这里
我们可以使用线性筛处理:
int pri[N],cnt;
int vis[N];
int sp1[N],sp2[N];//sp1存F_1的sp,sp2存F_2的sp
void init(int maxn)
{
for(int i=2;i<=maxn;i++)
{
if(!vis[i])
{
pri[++cnt] = i;
sp1[cnt] = (sp1[cnt-1] + i)%mod;
sp2[cnt] = (sp2[cnt-1] + i * i)%mod;
}
for(int j=1;j<=cnt;j++)
{
if(i * pri[j] > maxn)break;
vis[i * pri[j]] = 1;
if(i % pri[j] == 0)break;
}
}
}
接下来处理
由于我们使用
所以第二个参数可以滚动掉。
但是由于
我们发现,
只有
但是如果直接存,下标还是到
所以考虑这样的储存策略:
使用离散化的方法,使用序号
使用
当数字
所以令
因为存在
反之则
因为存在
代码如下:
其中
int f1(int x)
{
x%=mod;
return x * (x + 1)/2 % mod;
}
int f2(int x)
{
x%=mod;
return x * (x+1) % mod * ((2 * x + 1)% mod) % mod * inv6 % mod;
}
int w[N],ind1[N],ind2[N],tot,g1[N],g2[N];
for(int l = 1,r;l<=n;l = r + 1) //枚举x可能出现的值(整除分块)
{
r = n / (n / l);
//w数组的含义:id(第几个可能的x值)转实际值
w[++tot] = n/l; //表示第tot个值是n/l
//计算g(x,0)
//实际处理中,我们一般不含f(1)(这里f(1)=1),只要最后加回即可
g1[tot] = f1(w[tot])-1;
g2[tot] = f2(w[tot])-1;
//ind1[x]为x的id
//ind2[x]为n/x的id
if(w[tot] <= sqr) // 当x<=sqrt(n)时,存于ind1,否则存于ind2
ind1[w[tot]] = tot;
else
ind2[n/w[tot]] = tot;
}
然后,根据公式递推
int getid(int x) // x to id
{
if(x <= sqr)
return ind1[x];
else
return ind2[n/x];
}
for(int i=1;i<=cnt;i++)
{
for(int j=1;j<=tot && pri[i] * pri[i] <= w[j];j++)
{
g1[j] =(g1[j] - pri[i] * (g1[getid(w[j]/pri[i])] - sp1[i-1]) % mod + mod)%mod;
g2[j] =(g2[j] - pri[i] * pri[i] % mod * (g2[getid(w[j]/pri[i])] - sp2[i-1]) % mod+mod)%mod;
}
}
最后计算
int s(int x,int k)
{
if(pri[k] > x)
return 0;
int ans = ((g2[getid(x)] - g1[getid(x)] + mod)%mod - (sp2[k] - sp1[k] + mod)%mod + mod)%mod;
for(int i=k+1;i<=cnt && pri[i] * pri[i] <= x;i++)
{
for(int a = 1,p = pri[i];p <= x;p *= pri[i],a++)
{
ans =(ans + p % mod * ((p% mod-1) % mod * (s(x/p,i) + (a!=1)) % mod))%mod;
}
}
return ans;
}
汇总代码:
#include<bits/stdc++.h>
#define int long long
const int N = (int)2e5+7,mod = (int)1e9+7;
const int inv2 = 500000004,inv6 = 166666668;
using namespace std;
int pri[N],cnt;
int vis[N];
int sp1[N],sp2[N];
void init(int maxn)
{
for(int i=2;i<=maxn;i++)
{
if(!vis[i])
{
pri[++cnt] = i;
sp1[cnt] = (sp1[cnt-1] + i)%mod;
sp2[cnt] = (sp2[cnt-1] + i * i)%mod;
}
for(int j=1;j<=cnt;j++)
{
if(i * pri[j] > maxn)break;
vis[i * pri[j]] = 1;
if(i % pri[j] == 0)break;
}
}
}
int n,sqr;
int w[N],ind1[N],ind2[N],tot,g1[N],g2[N];
int getid(int x)
{
if(x <= sqr)
return ind1[x];
else
return ind2[n/x];
}
int f1(int x)
{
x%=mod;
return x * (x + 1)/2 % mod;
}
int f2(int x)
{
x%=mod;
return x * (x+1) % mod * ((2 * x + 1)% mod) % mod * inv6 % mod;
}
int s(int x,int k)
{
if(pri[k] > x)
return 0;
int ans = ((g2[getid(x)] - g1[getid(x)] + mod)%mod - (sp2[k] - sp1[k] + mod)%mod + mod)%mod;
for(int i=k+1;i<=cnt && pri[i] * pri[i] <= x;i++)
{
for(int a = 1,p = pri[i];p <= x;p *= pri[i],a++)
{
ans =(ans + p % mod * ((p% mod-1) % mod * (s(x/p,i) + (a!=1)) % mod))%mod;
}
}
return ans;
}
signed main()
{
cin>>n;
sqr = sqrt(n);
init(sqr);
for(int l = 1,r;l<=n;l = r + 1)
{
r = n / (n / l);
w[++tot] = n/l;
g1[tot] = f1(w[tot])-1;
g2[tot] = f2(w[tot])-1;
if(w[tot] <= sqr)
ind1[w[tot]] = tot;
else
ind2[n/w[tot]] = tot;
}
for(int i=1;i<=cnt;i++)
{
for(int j=1;j<=tot && pri[i] * pri[i] <= w[j];j++)
{
g1[j] =(g1[j] - pri[i] * (g1[getid(w[j]/pri[i])] - sp1[i-1]) % mod + mod)%mod;
g2[j] =(g2[j] - pri[i] * pri[i] % mod * (g2[getid(w[j]/pri[i])] - sp2[i-1]) % mod+mod)%mod;
}
}
//f(1)=1
cout<<(s(n,0)+1)%mod;
return 0;
}
于是,你就拿到了
对本例题的汇总
min25 筛的过程中,需要求解的有:
函数 函数 函数
对于
查看整个 min25 筛的过程,使用到的
所以预处理出
对于
对于
然后根据推导出的公式递推。
对于
0x21 补充例题
Loj6235 区间素数个数
首先,我们按照需求,构造出需要求得的函数。
如果你对 min25 筛的过程足够熟悉,你就会发现
的相似之处。
不妨变换一下:
于是我们很快就打出代码:
#include <bits/stdc++.h>
#define ll long long
using namespace std;
//变量名的意义和之前时一样的
ll n, sqr,w[1000005],g[1000005];
int pri[1000005], cnt;
int vis[1000005];
int ind1[1000005], ind2[1000005], tot;
void init(int maxn) {
//k值为0,sp[x]即为1~p_x的质数个数
for (int i = 2; i <= maxn; i++) {
if (!vis[i]) {
pri[++cnt] = i;
//sp[cnt] = cnt;
//没有必要多用一个数组储存sp
}
for (int j = 1; j <= cnt; j++) {
if (i * pri[j] > maxn)
break;
vis[i * pri[j]] = 1;
if (i % pri[j] == 0)
break;
}
}
}
int getid(ll x) {
if (x <= sqr)
return ind1[x];
else
return ind2[n / x];
}
signed main() {
scanf("%lld",&n);
sqr = sqrt(n);
init(sqr);
for (ll l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
w[++tot] = n / l;
g[tot] = w[tot] - 1;
if (w[tot] <= sqr) {
ind1[w[tot]] = tot;
} else {
ind2[n / w[tot]] = tot;
}
}
for (int j = 1; j <= cnt; j++) {
for (int i = 1; i <= tot; i++) {
//k=0
if ((ll)pri[j] * pri[j] > w[i])break;
g[i] -= (g[getid(w[i] / pri[j])] - (j-1));
}
}
printf("%lld",g[getid(n)]);
}
100pts
但请记住,在这里:
for (int j = 1; j <= cnt; j++) {
for (int i = 1; i <= tot; i++) {
//k=0
if ((ll)pri[j] * pri[j] > w[i])break;
g[i] -= (g[getid(w[i] / pri[j])] - (j-1));
}
}
请记住别写成这样(60pts TLE 的痛):
for (int j = 1; j <= cnt; j++) {
for (int i = 1; i <= tot; i++) {
//k=0
if ((ll)pri[j] * pri[j] <= w[i])
g[i] -= (g[getid(w[i] / pri[j])] - (j-1));
}
}