[学习笔记] Min_25筛
这个东西 \(\tt jzm\) 又会了！但是我自己看了一会还是看懂了。
什么是Min25？
据说是一个叫 \(\tt Min25\) 的人发明的，跟杜教筛的玄学程度有得一拼，不过并不需要很多的前置知识。
它是用来解决这样一类问题：已知 \(f(p^k)\) 是关于质数 \(p\) 的多项式，\(f(x)\) 是积性函数，求 \(\sum_{i=1}^n f(i)\) ，做法是把积性函数拆成若干个完全积性函数，算出对于 \(n\) 以内质数的函数值，再用递归的形式算答案，总的来说，体现了用 函数结构，找递归子问题 来加速计算的思想（再杜教筛的题里面也很常见）。
如果你以前没有学过 \(\tt Min25\) ，你肯定听不懂我在讲什么，没关系，我们慢慢来解析这个算法。
由于 \(1\) 很特殊，所以我们下面的算法是不考虑 \(1\) 的，最后把 \(1\) 的答案加上即可。
分类
首先对 \(i\) 按照质数和合数分类：
对于后面的数枚举 最小质因数 ，因为合数的最小质因数 \(\leq\sqrt n\) 所以会快（这也印证了为什么单独讨论质数）：
其中 \(minp\) 表示 \(i\) 的最小质因数，发现了吗？我们把质数提出来后，利用积性函数的性质，把最小质因子拆出来，就得到了一个子问题 ，这时候我们的思路就渐渐明晰了，算法分为两部分：处理前面质数的答案，处理后面的递归子问题。
第一部分
这是 \(\tt Min25\) 的神奇之处，我们可以在不知道所有质数的情况下求出所有质数的答案之和。
设 \(g(n,j)\) 表示 \(n\) 以内的数，满足是个质数或者最小质因数 \(>p_j\) 的合数的函数值，可能我们需要的函数值是一个多项式，但是这里我们把他改成一个形如 \(x^k\) 的东西，这样他就是一个 完全积性函数 ，因为 \(g(..,0)\) 是好算的（直接求和），然后我们用 \(\sqrt n\) 以内的质数慢慢 把合数筛掉 ，最后就只剩质数了，我们来看一看怎么筛的：
减掉的就是包含 \(p_j\) 这个质因子，且它是最小质因子的数，因为是完全积性函数所以不用把所有的质因子提出来，而只用提一个。后面减掉的是 \(p_{j-1}\) 以内的质数，这些是不应该算进去的，就是 \(sp\) 数组，它可以线性筛出来。
第一维很大怎么办，但是你发现我们用到的只有 \(\frac{n}{x}\) 这些取值，一共只有 \(\sqrt n\) 个，所以会很快。
第二部分
要解决递归的问题了，设 \(S(n,j)\) 表示 \(n\) 以内的数最小质因子大于 \(p_j\) 的函数值之和，假设我们已经做完了第一部分得到了 \(g(i)\) 表示 \(i\) 以内质数的函数值之和，那么就这样递归：
后面的 \([e\not=1]\) 就表示如果不是只提了一次，那么就是某个质数大于 \(1\) 的幂次，这部分的答案是我们没有统计到的，所以加上。这一部分的复杂度我不太清楚，所以跟杜教筛一样玄学啊。
例题
这道题就是把原函数值拆成 \(x^2,x^3\) 的完全积性函数 \(g\)，其他就是正常写法啦。
今天写 \(\tt Min25\) 的时候遇到了问题，一开始筛 \(g\) 的时候要注意只能这么写 p[i]*p[i]<=w[j] ，必须保证有合数可以筛要不然会出问题
#include <cstdio>
#include <iostream>
#include <cmath>
using namespace std;
const int N = 200000;
const int M = N+5;
const int MOD = 1e9+7;
const int inv = 166666668;
#define int long long
int read()
{
int x=0,f=1;char c;
while((c=getchar())<'0' || c>'9') {if(c=='-') f=-1;}
while(c>='0' && c<='9') {x=(x<<3)+(x<<1)+(c^48);c=getchar();}
return x*f;
}
int n,cnt,tot,p[M],vis[M],w[M],sp1[M],sp2[M];
int sqr,id1[M],id2[M],g1[M],g2[M];
void init(int n)
{
for(int i=2;i<=n;i++)
{
if(!vis[i])
{
p[++cnt]=i;
sp1[cnt]=(sp1[cnt-1]+i)%MOD;
sp2[cnt]=(sp2[cnt-1]+i*i)%MOD;
}
for(int j=1;j<=cnt && i*p[j]<=n;j++)
{
vis[i*p[j]]=1;
if(i%p[j]==0) break;
}
}
}
int S(int x,int y)
{
if(x<=p[y]) return 0;
int k=x<=sqr?id1[x]:id2[n/x];
int ans=(g2[k]-g1[k]-(sp2[y]-sp1[y]))%MOD;
for(int i=y+1;i<=cnt && p[i]*p[i]<=x;i++)
{
int pr=p[i];
for(int e=1;pr<=x;e++,pr*=p[i])
{
int xx=pr%MOD;
ans=(ans+xx*(xx-1)%MOD*(S(x/pr,i)+(e!=1)))%MOD;
}
}
return ans;
}
signed main()
{
init(N);
n=read();sqr=sqrt(n);
for(int l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
w[++tot]=n/l;
g1[tot]=w[tot]%MOD;//要模要模
g2[tot]=g1[tot]*(g1[tot]+1)%MOD*(2*g1[tot]+1)%MOD*inv%MOD-1;
g1[tot]=g1[tot]*(g1[tot]+1)/2%MOD-1;
if(n/l<=sqr) id1[n/l]=tot;//编号方式是x
else id2[n/(n/l)]=tot;//编号方式是n/x
}
for(int i=1;i<=cnt;i++)
for(int j=1;j<=tot && p[i]*p[i]<=w[j];j++)
{
int k=w[j]/p[i]<=sqr?id1[w[j]/p[i]]:id2[n/(w[j]/p[i])];
g1[j]=(g1[j]-p[i]*(g1[k]-sp1[i-1]))%MOD;
g2[j]=(g2[j]-p[i]*p[i]%MOD*(g2[k]-sp2[i-1]))%MOD;
}
printf("%lld\n",(S(n,0)+1+MOD)%MOD);
}
