既然选择了远方,便只顾风雨兼行|

H_W_Y

园龄:1年11个月粉丝:28关注:15

Number Theory(5)

12 奇妙筛法

已开工。


12.1 杜教筛

12.2 Min-25 筛

Min_25 筛可以在 O(n34logn) 的时间复杂度内解决一类 积性函数 前缀和问题。

说形象点就是 1s 能跑 1010,相当优秀。

要求:

  • f(p) 是关于 p 的低阶多项式,或者可以快速求出。
  • f(pc) 可以快速求出。

以下 p 代表质数,pi 表示第 i 小的质数,而 p0=1

同时我们钦定 lpf(i) 表示 i 的最小质因子。


首先,我们尝试对其进行推导,分成质数和合数两个部分

i=1nf(i)=1pnf(p)+1pen,1pnf(pe)(1inpe,lpf(i)>pf(i))

前面的部分是枚举素数,而后面的部分就是对于每一个合数枚举了其最小的质因子以及其次数。


考虑第一块的质数部分如何求值?

由于我们在前面要求了 f(p) 是关于 p 的低阶多项式,也就是

f(p)=iaipci

那么

1pnf(p)=iai1pnpci

现在问题就变成了求 1pnpk 了。

我们只对一个 k 拿出来讨论,如果多个就分别做以下的操作就可以了。


我们尝试用一个神秘的 dp,设

g(n,j)=i=1n[iPlpf(i)>pj]ik

(至于它有什么用,你先别急)

我们考虑 g(n,j) 是如何转移过来的。

不难想到用 g(n,j1) 进行转移,也就是减去了 lpf(i)=pj 的部分,它可以被表示成

pjkg(npj,j1)

由于 ik 是完全积性函数,我们不需要讨论前后是否互质!

而仔细一想,发现我们把质数多算了一次,它一共出现了 g(pj1,j1) 次,那么实际上减去的就是

pjk(g(npj,j1)g(pj1,j1))

其实后面的 g(pj1,j1) 就是前 j1 个质数的 k 次方和。

这样一来,我们就有了递推式

g(n,j)=g(n,j1)pjk(g(npj,j1)g(pj1,j1))

发现这是可以依次 dp 下去的,只要我们知道

i=2nik

就可以用质数处理出来。(但是时间复杂度根本过不去啊?你先别急)


回到引入 g 时的问题,它有什么用呢?

我们发现如果 px 是最大的 n 的质数,那么 g(n,x) 就是我们想要的答案。

所以如果通过 dp 推出 g(n,x) 的值就完全胜利,而我们把它简记成 g(n)。(具体细节你先别急)


假设我们现在已经可以快速算出 g 了,那原式会变成什么样子呢?

发现后面的那一块也是不好处理的,所以我们采用类似的 dp 方式,设 S(n,x) 表示 1nlpf(i)>px 的函数值和,容易发现答案就是 S(n,0),而我们有

S(n,x)=g(n)spx+pken,x<kf(pke)(S(npke,k)+[e>1])

其中 spx 表示前 x 个素数的答案,而容易发现当 px>n 是没有意义的,因为只会剩下一些大质数,所以我们的 pxn 级别的,完全可以预处理出来。

而后面的 [e>1] 则表示当前的 i=pe 的情况。

这样的递归是可以直接求的,时间复杂度为 O(n34logn) 的。具体证明可见 OI-wiki


现在我们回到对 g 的求解,具体怎么求?

发现直接处理出 1n 的所有值显然不太可能,n 太大了,但仔细观察两个式子,发现每次需要用到的就是

g(npj),g(npke)

这不就是类似于整除分块的东西么?

我们在第一章证明过

nab=nab

所以在整个过程中只会用到 1,2,,n,nn,,n 位置的值,就是整除分块的 ni 那些值。

在第九章我们知道个数是 2n 的,于是我们只需要对这 O(n)g 值进行 dp 预处理就可以了。


再来讲一讲具体的实现细节,我们用 id1[]id2[] 分别记录 ini>n 的值在 g 数组中的下标,而对于 i>n 的部分直接维护 ni 就可以了。

同时,由于 1 既不是质数也不是合数,所以在上面的推到中它都不应该被算进去,也就是 gS 都把 1 抠掉了,所以最后的答案是 S(n,0)+1


于是就可以通过 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;
}

现在我们换几个函数来筛,也就是我们最熟悉的 μφ

它们在质数处的取值也是相当容易的,μ(p)=1,φ(p)=p1,这样的低阶多项式是可以较为方便计算的。

直接把两个分开写出来就是这样的

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,最下面的是杜教筛)

image

于是我们开始了常数优化。


优化也是简单的,你把预处理放在一起,就可以跑到 2.95s:

image

进而,把两个函数一起做(再把不需要 long long 的地方换成 int),于是就可以做到 1.44s!!!

image

代码

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 阶与原根

n(n>1) 的简化剩余系在模 n 乘法意义下封闭且构成群。容易证明其满足封闭性和结合律,存在逆元和单位元。将群相关定义应用在其上,得到阶和原根(循环群和生成元)的概念。

                                 ——Alex_Wei

13.1 阶

由欧拉定理可知,对于 aZ,mN,若 am,则有

aφ(m)1(modm)

因此满足同余式 an1(modm) 的最小正整数 n 存在,这个 n 称为 am 的阶,记作 δm(a)ordm(a)


阶有一些性质。

a,a2,,aδm(a)m 两两不同余。

这和我们在 3.5 里面探讨的遗留问题较为相似。

考虑反证,假设存在 ij,且 aiaj(modm),则有 a|ij|1(modp)

但是显然有 0<|ij|<δm(a),与阶的定义矛盾。得证。


an1(modm),则 δm(a)n

n=δm(a)q+r,0r<δm(a)

r>0,则

arar(aδm(a))qan1(modm)

与阶的定义矛盾(最小性)。得证。


mN,a,bZ,am,bm,则

δm(ab)=δm(a)δm(b)

的充分必要条件是

δm(a)δm(b)

  • 必要性:因为 aδm(a)1(modm)bδm(b)1 可知

    (ab)lcm(δm(a),δm(b))1(modm)

    由于上面的性质,我们知道

    δm(ab)lcm(δm(a),δm(b))

    所以我们有

    δm(a)δm(b)lcm(δm(a),δm(b))

    gcd(δm(a),δm(b))=1

  • 充分性:因为

    (ab)δm(ab)1(modm)

    那么

    1(ab)δm(ab)δm(b)aδm(ab)δm(b)(modm)

    所以

    δm(a)δm(ab)δm(b)

    又由于 gcd(δm(a),δm(b))=1

    δm(a)δm(ab)

    同理可得

    δm(b)δm(ab)

    所以

    δm(a)δm(b)δm(ab)

    而我们有

    (ab)δm(a)δm(b)(aδm(a))δm(b)×(bδm(b))δm(a)1(modm)

    δm(ab)δm(a)δm(b)

    所以

    δm(ab)=δm(a)δm(b)

得证。


kN,mN,aZ,gcd(a,m)=1,则

δm(ak)=δm(a)gcd(δm(a),k)

很重要的一个性质。

注意到

akδm(ak)=(ak)δm(ak)1(modm)δm(a)kδm(ak)δm(a)gcd(δm(a),k)δm(ak)

而由于 aδm(a)1(modm)

(ak)δm(a)gcd(δm(a),k)=(aδm(a))kgcd(δm(a),k)1(modm)

所以

δm(ak)δm(a)gcd(δm(a),k)

δm(ak)=δm(a)gcd(δm(a),k)

得证。


容易发现其实上面的证明方法都挺类似的。/kk


那如何求 δm(a) 呢?(其中 ma

  • 直接使用 BSGS 求 ax1(modm) 的最小正整数解。时间复杂度 O(m)

  • 根据阶的性质,对于 x=kδm(a),必然有 ax1(modm)

    因此,我们考虑先让 x=φ(m),然后对于 φ(m) 的每个质因子 p,用 x 不断试除 p 知道 x 无法整除 p 或者 aap1(modm)

    时间复杂度 O(m+log2m),因为需要分解质因数和求解欧拉函数,用 Pollard_rho 可以优化为 m14

    m 不大,则可以 O(m) 线性处理每个数最小质因子,单次查询即可 log2m


13.2 原根

原根——数论的神!


mN,gZ。若 gcd(g,m)=1,且 δm(g)=φ(m),则称 g 为模 m 的原根。

我们需要原根的原因是它可以 生成m 的缩系,即它的若干次幂取遍了与 m 互质的数。这使得我们在考虑解与模数互质且模数存在原根的同余方程时可以用原根的若干次幂代替未知数。这一点在求解高次剩余时非常有用。


原根判定定理:设 m3,gcd(g,m)=1,则 g 是模 m 的原根的充要条件是,对于 φ(m) 的每个素因子 p,都有 gφ(m)p1(modm)

这和我们上面讲的求阶的方法非常像。

  • 必要性:显然。

  • 充分性:因为 φ(m)p 的所有因子取遍了 φ(m) 除了其本身的所有因子,所以若存在 ad1(modm),使得 d<φ(m),则必然存在 φ(m) 的质因子 p 使得 dφ(m)p,这就说明存在 aφ(m)p1(modm),与假设矛盾。得证


原根个数:若一个数 m 有原根,则它原根的个数为 φ(φ(m))

m 有原根 g,则:

δm(gk)=δm(g)gcd(δm(g),k)=φ(m)gcd(φ(m),k)

所以如果 gcd(φ(m),k)=1,则有 δm(gk)=φ(m),也就是 gk 也是模 m 的原根。

而满足 gcd(φ(m),k)=11kφ(m)kφ(φ(m)) 个。得证。


原根存在定理:一个数存在原根当且仅当 m=2,4,pa,2pa,其中 p 是奇素数,aN

可以当一个结论背住,具体证明比较复杂,可见 OI-wiki


最小原根的范围估计

在 1959 年,中国数学家王元证明了最小原根的级别为 O(m0.25+ϵ)。结合原根存在定理和原根判定定理,我们得以在 O(m) 预处理每个数的最小质因子后单次 O(m0.25ω(m)logm) 求出 m 的最小原根,从而求出 m 的所有原根。


关于原根,还有一个有趣且重要的性质:

m 存在原根 g,则 n=φ(m),则 m 的简化剩余系与 n 次单位根同构。

这是因为对于任何一个 am,其均可以写成 gk 的形式。而任何一个 n 次单位根 ωk 也可以写成 ωnk 的形式。

同时,gk1(modp),且 ωnn=1。因此,我们可以将单位圆应用在 m 的简化剩余系上。


这也是为什么对于特定模数,我们可以使用原根代替 n 次单位根进行 FFT,所以当需要使用 2k 次单位根时,若 k 不超过 φ(p) 当中质因子 2 的个数且 p 存在原根,即可用原根代替 2k 次单位根。

于是就得到了 NTT。


13.3 例题

放两道简单题。

P5605 小 A 与两位神仙

这个游戏的规则是这样的:两个神仙先选定一个正整数 m,保证 m 是一个奇质数 p 的正整数次幂。然后进行 n 轮游戏,每轮中造梦者选定一个正整数 x,杰瑞米选定一个正整数 y,保证 (x,m)=1,(y,m)=1,即 xm 互质,ym 互质,接下来询问小 A 是否存在非负整数 a 使得 xay(modm)

神仙们说小 A 只有在每一轮游戏中都回答正确才能回到正常的生活中,不得已之下他只好求助于聪明的你。

1n2×1043m10181x,y<m

注意到模数 p 是奇素数的幂,所以存在原根 g


我们尝试将 xy 表示成 gXgY 的形式。


14 高次剩余问题

15 卢卡斯定理



后记

参考文献:

OI-Note Chapter4.1 整除与同余

OI-Note Chapter4.2 数论函数与求和

初等数论学习笔记 I:同余相关

初等数论学习笔记 II:分解质因数

初等数论学习笔记 III:数论函数与筛法

Pollard_Rho 算法 - cjTQX - 博客园 (cnblogs.com)

《具体数学》第四章

本文作者:H_W_Y

本文链接:https://www.cnblogs.com/H-W-Y/p/18203768/NumberTheory5

版权声明:本作品采用知识共享署名-非商业性使用-禁止演绎 2.5 中国大陆许可协议进行许可。

posted @   H_W_Y  阅读(64)  评论(0编辑  收藏  举报
点击右上角即可分享
微信分享提示
评论
收藏
关注
推荐
深色
回顶
收起