Number Theory(5)
12 奇妙筛法
已开工。
12.1 杜教筛
12.2 Min-25 筛
Min_25 筛可以在
说形象点就是
要求:
是关于 的低阶多项式,或者可以快速求出。 可以快速求出。
以下
同时我们钦定
首先,我们尝试对其进行推导,分成质数和合数两个部分
前面的部分是枚举素数,而后面的部分就是对于每一个合数枚举了其最小的质因子以及其次数。
考虑第一块的质数部分如何求值?
由于我们在前面要求了
那么
现在问题就变成了求
我们只对一个
我们尝试用一个神秘的 dp,设
(至于它有什么用,你先别急)
我们考虑
不难想到用
由于
而仔细一想,发现我们把质数多算了一次,它一共出现了
其实后面的
这样一来,我们就有了递推式
发现这是可以依次 dp 下去的,只要我们知道
就可以用质数处理出来。(但是时间复杂度根本过不去啊?你先别急)
回到引入
我们发现如果
所以如果通过 dp 推出
假设我们现在已经可以快速算出
发现后面的那一块也是不好处理的,所以我们采用类似的 dp 方式,设
其中
而后面的
这样的递归是可以直接求的,时间复杂度为
现在我们回到对
发现直接处理出
这不就是类似于整除分块的东西么?
我们在第一章证明过
所以在整个过程中只会用到
在第九章我们知道个数是
再来讲一讲具体的实现细节,我们用
同时,由于
于是就可以通过 P5325 【模板】Min_25 筛 了。代码。
const ll H=1e9+7;
int cnt=0,tot=0;
ll p[N],g1[N],g2[N],s1[N],s2[N],w[N],id1[N],id2[N],B,n,i6,i2;
bool vis[N];
ll qpow(ll a,ll b){
ll res=1;
while(b){if(b&1) res=res*a%H;a=a*a%H,b>>=1;}
return res;
}
ll S1(ll x){x%=H;return x*(x+1)%H*i2%H;}
ll S2(ll x){x%=H;return x*(x+1)%H*(2*x+1)%H*i6%H;}
ll pos(ll x){return x<=B?id1[x]:id2[n/x];}
ll S(ll x,int y){
if(x<=p[y]) return 0;
int k=pos(x),res=((g2[k]-g1[k]-s2[y]+s1[y])%H+H)%H;
for(int i=y+1;i<=cnt&&p[i]*p[i]<=x;i++)
for(ll j=1,nw=p[i];nw<=x;++j,nw*=p[i]){
ll z=nw%H;
res=(res+(z-1)*z%H*(S(x/nw,i)+(j>1))%H)%H;
}
return res;
}
int main(){
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
cin>>n,B=sqrt(n);
i2=qpow(2,H-2),i6=qpow(6,H-2);
for(int i=2;i<=B;i++){
if(!vis[i]){
p[++cnt]=i;
s1[cnt]=(s1[cnt-1]+i)%H;
s2[cnt]=(s2[cnt-1]+1ll*i*i%H)%H;
}
for(int j=1;j<=cnt&&i*p[j]<=B;j++){
vis[i*p[j]]=true;
if(i%p[j]==0) break;
}
}
for(ll l=1,r;l<=n;l=r+1){
r=min(n/(n/l),n);
w[++tot]=n/l;
g1[tot]=S1(n/l)-1,g2[tot]=S2(n/l)-1;
if(n/l<=B) id1[n/l]=tot;
else id2[n/(n/l)]=tot;
}
for(int i=1;i<=cnt;i++){//g 数组的 dp,用了滚动数组
for(int j=1;p[i]*p[i]<=w[j]&&j<=tot;j++){
int cur=pos(w[j]/p[i]);
g1[j]=(g1[j]-p[i]*(g1[cur]-s1[i-1])%H+H)%H;
g2[j]=(g2[j]-p[i]*p[i]%H*(g2[cur]-s2[i-1])%H+H)%H;
}
}
cout<<(S(n,0)+1)%H;
return 0;
}
现在我们换几个函数来筛,也就是我们最熟悉的
它们在质数处的取值也是相当容易的,
直接把两个分开写出来就是这样的
const int N=1e6+5;
int cnt=0,tot=0;
ll p[N],w[N],id1[N],id2[N],n,B;
bool vis[N>>2];
int pos(ll x){return x<=B?id1[x]:id2[n/x];}
void Pr(int n){
memset(vis,0,sizeof(vis));
for(int i=2;i<=n;i++){
if(!vis[i]) p[++cnt]=i;
for(int j=1;j<=cnt&&i*p[j]<=n;j++){
vis[i*p[j]]=true;
if(i%p[j]==0) break;
}
}
}
namespace Mu{
ll g[N];
ll S(ll x,int y){
if(x<=p[y]) return 0;
ll res=-g[pos(x)]+y;
for(int i=y+1;i<=cnt&&p[i]*p[i]<=x;i++)
res-=S(x/p[i],i);
return res;
}
void SOL(){
for(int i=1;i<=tot;i++) g[i]=w[i]-1;
for(int i=1;i<=cnt;i++){
for(int j=1;p[i]*p[i]<=w[j]&&j<=tot;j++){
int cur=pos(w[j]/p[i]);
g[j]+=i-g[cur]-1;
}
}
cout<<(S(n,0)+1)<<'\n';
}
}
namespace Phi{
ll g1[N],g2[N],s[N];
ll S(ll x,int y){
if(x<=p[y]) return 0;
int k=pos(x);
ll res=g1[k]-g2[k]-(s[y]-y);
for(int i=y+1;i<=cnt&&p[i]*p[i]<=x;i++){
for(ll j=1,nw=1;nw*p[i]<=x;++j,nw*=p[i])
res+=nw*(p[i]-1)*(S(x/(nw*p[i]),i)+(j>1));
}
return res;
}
void SOL(){
for(int i=1;i<=cnt;i++) s[i]=s[i-1]+p[i];
for(int i=1;i<=tot;i++) g1[i]=1ll*w[i]*(w[i]+1)/2-1,g2[i]=w[i]-1;
for(int i=1;i<=cnt;i++)
for(int j=1;p[i]*p[i]<=w[j]&&j<=tot;j++){
int cur=pos(w[j]/p[i]);
g1[j]-=p[i]*(g1[cur]-s[i-1]);
g2[j]+=i-1-g2[cur];
}
cout<<(S(n,0)+1)<<' ';
}
}
void SOLVE(){
cin>>n,B=sqrt(n),cnt=tot=0;
memset(vis,0,sizeof(vis));
Pr(B);
for(ll l=1,r;l<=n;l=r+1){
r=min(n,n/(n/l));
w[++tot]=n/l;
if(n/l<=B) id1[n/l]=tot;
else id2[n/(n/l)]=tot;
}
Phi::SOL(),Mu::SOL();
}
int main(){
/*2024.4.26 H_W_Y P4213 【模板】杜教筛 Min_25 筛*/
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
int _;cin>>_;
while(_--) SOLVE();
return 0;
}
但结果并不理想:(前两个是 Min_25,最下面的是杜教筛)
于是我们开始了常数优化。
优化也是简单的,你把预处理放在一起,就可以跑到 2.95s:
进而,把两个函数一起做(再把不需要 long long 的地方换成 int),于是就可以做到 1.44s!!!
代码。
const int N=1e6+5;
int cnt=0,tot=0,p[N],n,B,w[N],id1[N],id2[N],s[N];
ll g1[N],g2[N];
bool vis[N>>2];
inline int pos(int x){return x<=B?id1[x]:id2[n/x];}
inline void Pr(int n){
for(int i=2;i<=n;i++){
if(!vis[i]) p[++cnt]=i,s[cnt]=s[cnt-1]+i;
for(int j=1;j<=cnt&&i*p[j]<=n;j++){
vis[i*p[j]]=true;
if(i%p[j]==0) break;
}
}
}
inline pll S(int x,int y){
if(x<=p[y]) return {0,0};
int k=pos(x);
pll res={g1[k]-g2[k]-(s[y]-y),-g2[k]+y};
for(int i=y+1;i<=cnt&&p[i]*p[i]<=x;i++)
for(int j=1,nw=1;1ll*nw*p[i]<=x;++j,nw*=p[i]){
pll cur=S(x/(nw*p[i]),i);
if(nw==1) res.se-=cur.se;
res.fi+=1ll*nw*(p[i]-1)*(cur.fi+(j>1));
}
return res;
}
inline void SOLVE(){
cin>>n,B=sqrt(n),tot=0;
for(ll l=1,r;l<=n;l=r+1){
r=n/(n/l);
w[++tot]=n/l;
g1[tot]=1ll*(n/l)*(n/l+1)/2-1,g2[tot]=n/l-1;
if(n/l<=B) id1[n/l]=tot;
else id2[n/(n/l)]=tot;
}
for(int i=1;i<=cnt&&p[i]<=B;i++)
for(int j=1;p[i]*p[i]<=w[j]&&j<=tot;j++){
int cur=pos(w[j]/p[i]);
g1[j]-=1ll*p[i]*(g1[cur]-s[i-1]);
g2[j]+=i-1-g2[cur];
}
pll res=S(n,0);
cout<<(res.fi+1)<<' '<<(res.se+1)<<'\n';
}
int main(){
/*2024.4.26 H_W_Y P4213 【模板】杜教筛 Min_25 筛*/
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
int _;cin>>_;Pr(46340);
while(_--) SOLVE();
return 0;
}
12.3 Powerful Number 筛
12.4 扩展埃氏筛
13 阶与原根
模
的简化剩余系在模 乘法意义下封闭且构成群。容易证明其满足封闭性和结合律,存在逆元和单位元。将群相关定义应用在其上,得到阶和原根(循环群和生成元)的概念。
——Alex_Wei
13.1 阶
由欧拉定理可知,对于
因此满足同余式
阶有一些性质。
模 两两不同余。
这和我们在 3.5 里面探讨的遗留问题较为相似。
考虑反证,假设存在
但是显然有
若
,则 。
设
若
与阶的定义矛盾(最小性)。得证。
设
,则 的充分必要条件是
-
必要性:因为
和 可知由于上面的性质,我们知道
所以我们有
即
。 -
充分性:因为
那么
所以
又由于
,同理可得
所以
而我们有
故
所以
得证。
设
,则
很重要的一个性质。
注意到
而由于
所以
故
得证。
容易发现其实上面的证明方法都挺类似的。/kk
那如何求
-
直接使用 BSGS 求
的最小正整数解。时间复杂度 。 -
根据阶的性质,对于
,必然有 。因此,我们考虑先让
,然后对于 的每个质因子 ,用 不断试除 知道 无法整除 或者 。时间复杂度
,因为需要分解质因数和求解欧拉函数,用 Pollard_rho 可以优化为 。若
不大,则可以 线性处理每个数最小质因子,单次查询即可 。
13.2 原根
原根——数论的神!
设
我们需要原根的原因是它可以 生成 模
原根判定定理:设
这和我们上面讲的求阶的方法非常像。
-
必要性:显然。
-
充分性:因为
的所有因子取遍了 除了其本身的所有因子,所以若存在 ,使得 ,则必然存在 的质因子 使得 ,这就说明存在 ,与假设矛盾。得证
原根个数:若一个数
若
所以如果
而满足
原根存在定理:一个数存在原根当且仅当
可以当一个结论背住,具体证明比较复杂,可见 OI-wiki。
最小原根的范围估计:
在 1959 年,中国数学家王元证明了最小原根的级别为
关于原根,还有一个有趣且重要的性质:
若
这是因为对于任何一个
同时,
这也是为什么对于特定模数,我们可以使用原根代替
于是就得到了 NTT。
13.3 例题
放两道简单题。
P5605 小 A 与两位神仙
这个游戏的规则是这样的:两个神仙先选定一个正整数
,保证 是一个奇质数 的正整数次幂。然后进行 轮游戏,每轮中造梦者选定一个正整数 ,杰瑞米选定一个正整数 ,保证 ,即 与 互质, 与 互质,接下来询问小 A 是否存在非负整数 使得 。 神仙们说小 A 只有在每一轮游戏中都回答正确才能回到正常的生活中,不得已之下他只好求助于聪明的你。
, , 。
注意到模数
我们尝试将
14 高次剩余问题
15 卢卡斯定理
后记
参考文献:
Pollard_Rho 算法 - cjTQX - 博客园 (cnblogs.com)
《具体数学》第四章
本文作者:H_W_Y
本文链接:https://www.cnblogs.com/H-W-Y/p/18203768/NumberTheory5
版权声明:本作品采用知识共享署名-非商业性使用-禁止演绎 2.5 中国大陆许可协议进行许可。
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步