标签:ret str 题目 pre printf 因子 部分 概念 证明
%%%gmh差不多一年前学会min_25筛
%%%某初一大佬似乎已经会了min_25筛
菜哭了
以下记\(P\)为素数集合,\(P(n)\)为所有小于等于\(n\)的素数的集合。
\(minp(x)\)表示\(x\)的最小质因子
这种什么筛之类的,多是求积性函数的前缀和的算法。
min_25筛能做的积性函数要满足这两个条件:
min_25筛具体分为两个部分:求素数的答案,和求所有数的答案(不包括\(1\)的答案,\(1\)的答案最后自己加上)。
先说素数部分怎么求。
分别枚举每一项的系数,不同的系数分开来算。
以下记素数\(x\)的函数为\(f(x)\)(指某一特定的系数下)
设\(g(n,j)=\sum_{i \in P(n) 或 minp(i)>p_j}f(i)\)
考虑从\(g(n,j-1)\)得到\(g(n,j)\)
想想\(g(n,j-1)\)得到\(g(n,j)\)要去掉些什么,根据概念可知要去掉所有满足\(minp(x)=p_j\)的合数\(x\)。
这个合数最小是\(p_j^2\),于是当\(n<p_j^2\)时,有\(g(n,j)=g(n,j-1)\)
古爷:这就是min_25筛的精髓!
有了这个东西,意味着大于\(\sqrt n\)的质数都没必要算,
\(j\)最多取到\(|P(\sqrt n)|\)就可以了。
当\(n\ge p_j^2\)时,考虑将\(p_j\)的倍数筛去。
由于是积性函数,所以可以将\(f(p_j)\)提出来,剩下的就是\(g(\lfloor\frac{n}{p_j} \rfloor,j-1)\)。
真的是吗?其实这样还算多了小于\(p_j\)的所有质数的贡献,即\(g(p_j-1,j-1)\),减去即可。
于是这个时候\(g(n,j)=g(n,j-1)+f(p_j)(g(\lfloor\frac{n}{p_j} \rfloor,j-1)-g(p_j-1,j-1))\)
注意一点,上面说“由于是积性函数”,这里应该是广义的积性函数(即便不是题目想要求的),也就是\(f(x)\),形式为\(f(x)=x^k\)
这个步骤下我们最终只是需要素数的贡献和,至于合数的贡献,我们把它当做广义积性函数来算。
合数的贡献只是中间值,计算到最后是会消掉的。
临界:当\(j=0\)的时候,直接用自然数幂和来计算。注意减去\(1\)的贡献。
显然,\(g(n,|P(n)|)\)就是所有小于等于\(n\)的质数的答案。
形象化地理解这个东西:
当从\(P(n,j-1)\)转移到\(P(n,j)\)的时候,将\(minp(x)=p_j\)的合数\(x\)去掉。
这个东西看起来很像埃氏筛法。
接下来是计算所有数的答案。
这里的\(F(x)\)直接表示\(x\)的贡献(不需要分具体是哪一位)。
设\(S(n,j)=\sum_{minp(i)\ge p_j} f(i)\)
先列出式子:
\(S(n,j)=g(n,|P(n)|)-\sum_{k<j} F(k)+\sum_{k\ge j,p_k^2\le n}\sum_{e>0,p_k^{e+1}\le n}(F(p_k^e)S(\lfloor \frac{n}{p_k^e}\rfloor,k+1)+F(p_k^{e+1}))\)
前面的\(g(n,|P(n)|)-\sum_{k<j} f(k)\)是所有大于等于\(p_j\)的素数的答案。
后面就是枚举最小质因子\(p_k\)及其系数\(e\),提出来之后计算。
\(S(\lfloor \frac{n}{p_k^e}\rfloor,k+1)\)如果不为\(0\),则至少要满足\(p_k\le\lfloor\frac{n}{p_k^e}\rfloor\),移项一下就得到了\(p_k^{e+1}\le n\)
最后面之所以要加上那个东西,是因为\(S(\lfloor \frac{n}{p_k^e}\rfloor,k+1)\)中并不包括\(p_k\),所以要单独计算\(p_k^{e+1}\)的贡献。
注意,这里的积性函数是“真”的积性函数
对于计算到的每个数,每个质因子都只会枚举一遍。
并不一定是广义积性函数。
答案就是\(S(n,1)+1\),后面那个是\(1\)的贡献。
观察前面的那个式子,感觉好丑啊。
为什么还要枚举\(k\)?直接枚举\(p_j\)的系数就好了啊!
\(S(n,j)=S(n,j+1)+F(p_j)+\sum_{e>0,p_j^{e+1}\le n} (F(p_j^e)S(\lfloor\frac{n}{p_j^e}\rfloor)+F(p_j^{e+1}))\)
哇哈哈看起来似乎连\(g\)都不用求
然后得到了古爷的指导:这个东西,\(O(n)\)起步。
后面的那一段求和只有\(p_j^2\le n\)的时候是要做的,但是前面的那个要做\(|P(n)|\)遍。
于是就有个简单的优化:
当\(p_j^2>n\)的时候,直接计算\(\sum_{p_k^2>n} F(p_k)\),加上之后退出。
这个东西用\(g\)来求。所以还是要求\(g\)啊……
待填坑。
观察一下式子,可以发现,要计算到的\(n\)都满足\(n\in\{\lfloor\frac{N}{i}\rfloor|i\in[1,N],i \in Z\}\)
\(N\)是题目给出的\(N\)。
自己证
可以发现这个集合的大小不超过\(2\sqrt N\)
考虑预处理出每个\(g(n,|P(n)|)\)。
简写为\(g_n\)
先考虑怎么存储。对于每个\(n\),可以按小于等于或者大于\(\sqrt N\)讨论,后者借助\(\lfloor\frac{N}{n}\rfloor\)来存储。可以证明\(\lfloor\frac{N}{n}\rfloor\)对于不同的\(n\)不重复。
接下来考虑模拟筛法来求:
枚举质数\(p_j\),然后计算当前这一层的\(g_n\)
具体见代码
求\(S\)的时候,不用记忆化,不用记忆化,不用记忆化!
另外,用前面那个方法就可以了,中间那个方法听说会慢,我自己脑补出来的那个方法和前面那个差不多(无法理解!为什么看起来优美那么多却没有地变快?)
例题见洛谷板题。
大众做法
using namespace std;
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <cassert>
#define M 100010
#define mo 1000000007
#define ll long long
ll qpow(ll x,ll y=mo-2){
ll r=1;
for (;y;y>>=1,x=x*x%mo)
if (y&1)
r=r*x%mo;
return r;
}
ll _n,sq;
int p[M],np,sf[M],sg1[M],sg2[M];
bool inp[M];
int lsp[M];
#define f(p) ((ll)(p)*((p)-1))
ll w[M*2],cnt;
int id1[M],id2[M];
#define id(d) (d<=sq?id1[d]:id2[_n/d])
ll g1[M*2],g2[M*2];
void init(ll n){
ll sq=sqrt(n);
ll inv2=qpow(2),inv6=qpow(6);
for (ll i=1;i<=n;++i){
ll j=n/i;
w[++cnt]=j;
(j<=sq?id1[j]:id2[n/j])=cnt;
i=n/j;
j%=mo;
g1[cnt]=(j*(j+1)%mo*inv2-1+mo)%mo;
g2[cnt]=(j*(j+1)%mo*(2*j+1)%mo*inv6-1+mo)%mo;
}
for (int j=1;j<=np;++j)
for (int i=1;i<=cnt && (ll)p[j]*p[j]<=w[i];++i){
ll d=w[i]/p[j],k=id(d);
g1[i]=((g1[i]-(ll)p[j]*(g1[k]-sg1[j-1]))%mo+mo)%mo;
g2[i]=((g2[i]-(ll)p[j]*p[j]%mo*(g2[k]-sg2[j-1]))%mo+mo)%mo;
}
}
ll S(ll n,int j){
ll tmp=id(n),r=(g2[tmp]-g1[tmp])-(j?sf[j-1]:0);
r%=mo,r+=r<0?mo:0;
if (j>sq)
return r;
for (int k=j;k<=np && (ll)p[k]*p[k]<=n;++k)
for (ll w=p[k];w*p[k]<=n;w*=p[k])
(r+=f(w%mo)%mo*S(n/w,k+1)+f(w*p[k]%mo))%=mo;
return r;
}
int main(){
scanf("%lld",&_n);
sq=sqrt(_n);
for (int i=2;i<=sq;++i){
if (!inp[i]){
p[++np]=i;
sf[np]=(sf[np-1]+f(i))%mo;
sg1[np]=(sg1[np-1]+i)%mo;
sg2[np]=(sg2[np-1]+(ll)i*i)%mo;
}
lsp[i]=np;
for (int j=1;j<=np && (ll)i*p[j]<=sq;++j){
inp[i*p[j]]=1;
if (i%p[j]==0)
break;
}
}
init(_n);
printf("%lld\n",(1+S(_n,1))%mo);
return 0;
}
我的辣鸡做法
using namespace std;
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <cassert>
#define M 100010
#define mo 1000000007
#define ll long long
ll qpow(ll x,int y=mo-2){
ll r=1;
for (;y;y>>=1,x=x*x%mo)
if (y&1)
r=r*x%mo;
return r;
}
ll _n,sq;
int p[M],np,sf[M],sg1[M],sg2[M];
bool inp[M];
int lsp[M];
#define f(p) ((ll)(p)*((p)-1))
ll w[M*2],cnt;
int id1[M],id2[M];
#define id(d) (d<=sq?id1[d]:id2[_n/d])
ll g1[M*2],g2[M*2];
void init(ll n){
ll sq=sqrt(n);
ll inv2=qpow(2),inv6=qpow(6);
for (ll i=1;i<=n;++i){
ll j=n/i;
w[++cnt]=j;
(j<=sq?id1[j]:id2[n/j])=cnt;
i=n/j;
j%=mo;
g1[cnt]=(j*(j+1)%mo*inv2-1+mo)%mo;
g2[cnt]=(j*(j+1)%mo*(2*j+1)%mo*inv6-1+mo)%mo;
}
for (int j=1;j<=np;++j)
for (int i=1;i<=cnt && (ll)p[j]*p[j]<=w[i];++i){
ll d=w[i]/p[j],k=id(d);
g1[i]=((g1[i]-(ll)p[j]*(g1[k]-sg1[j-1]))%mo+mo)%mo;
g2[i]=((g2[i]-(ll)p[j]*p[j]%mo*(g2[k]-sg2[j-1]))%mo+mo)%mo;
}
}
ll S(ll n,int j){
ll r=0;
if (j>np || (ll)p[j]*p[j]>n){
int l=min(np,j-1),t=id(n);
r+=g2[t]-g1[t]-sf[l];
return (r%mo+mo)%mo;
}
r=(f(p[j]%mo)+S(n,j+1))%mo;
ll w=p[j];
for (int e=1;w*p[j]<=n;++e,w*=p[j])
(r+=f(w%mo)%mo*S(n/w,j+1)+f(w*p[j]%mo))%=mo;
return r;
}
int main(){
freopen("in.txt","r",stdin);
scanf("%lld",&_n);
sq=sqrt(_n);
for (int i=2;i<=sq;++i){
if (!inp[i]){
p[++np]=i;
sf[np]=(sf[np-1]+f(i))%mo;
sg1[np]=(sg1[np-1]+i)%mo;
sg2[np]=(sg2[np-1]+(ll)i*i)%mo;
}
lsp[i]=np;
for (int j=1;j<=np && (ll)i*p[j]<=sq;++j){
inp[i*p[j]]=1;
if (i%p[j]==0)
break;
}
}
init(_n);
printf("%lld\n",(1+S(_n,1))%mo);
return 0;
}
标签:ret str 题目 pre printf 因子 部分 概念 证明
原文地址:https://www.cnblogs.com/jz-597/p/12891605.html