数论算法

数论算法

Reaepita Lv2
  2022-09-09 19:47:59 2022-09-09 19:47:59 Created   2023-08-06 15:21:37 2023-08-06 15:21:37 Updated    

数论算法

快速乘

1
2
3
4
5
ll mul(ll a,ll b,ll mod)
{
	ll tmp=a*b-(ll)((long double)a/mod*b+0.5)*mod;
	return tmp<0?tmp+mod:tmp;
}

同余系相关

Exgcd

存在 因为所以

所以

递归求解，终止时必然有 。

解集为 

1
2
3
4
5
6
void exgcd(ll a,ll b,ll &x,ll &y)
{
    if(b==0){x=1,y=0;return;}
    exgcd(b,a%b,x,y);
    ll t=y;y=x-y*(a/b),x=t;
}

ExCRT

$$
\left{\right.
$$

每次合并两个方程 和

存在

求解二元一次方程

设

可得新的方程为

1
2
3
4
5
6
7
8
9
10
11
12
13
14
ll ExCRT(ll *a,ll *p,const int &n)
{
	ll A=a[1],M=p[1],k1,k2;
	for(int i=2;i<=n;i++)
	{
		ll c=((a[i]-A)%p[i]+p[i])%p[i],g=gcd(M,p[i]);
		if(c%g)return -1;
		exgcd(M,p[i],k1,k2);
		k1=mul(k1,c/g,p[i]);
		ll now=M/g*p[i];
		A=(mul(k1,M,now)+A)%now,M=now;
	}
	return (A%M+M)%M;
}

ExBSGS

根据欧拉定理 ， 与 互质。

可以在 的时间复杂度解决这个问题。

但是对于一般情况 不一定与 互质。

设

若 ，则存在

否则，若 解为 ， 无解

由 与 互质可得

继续分解直到 为止

这时候需要特判 是否存在与 同余的数。

因为 ，所以原问题就被转换为了 ，且 互质。

考虑分块求解这个问题，设

xxxxxxxxxx Point get_nearest_point_on_segment(Point a,Line a1){ Point b=a1.s,c=a1.e; double t=dot(a-b,c-b)/dot(c-b,c-b); if(dcmp(t)!=-1&&dcmp(1-t)!=-1)return Vector(b.x+(c.x-b.x)*t,b.y+(c.y-b.y)*t); if(dist(a,b)<dist(a,c))return b; return c;}cpp

只需要将所有的 存入Hash表中，再对所有 查询合法的值。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
int ExBSGS(int a,int b,int p)
{
	if(a==0&&b==0)return 1;
	if(a==0)return -1;
	if(b==1)return 0;
	unordered_map<int,int>vis;
	int g=gcd(a,p),mag=1,t=0;
	while(g!=1)
	{
		if(b%g)return -1;
		++t,p/=g,b/=g,mag=mag*1LL*(a/g)%p;
		if(mag==b)return t;
		g=gcd(a,p);
	}
	vis.clear();
	int n=sqrt(p)+1;
	for(int i=0,m=b;i<n;i++,m=m*1LL*a%p)vis[m]=i;
	a=Pow(a,n,p);
	for(int i=1,m=mag*1LL*a%p;i<=n;i++,m=m*1LL*a%p)if(vis.count(m))return i*n-vis[m]+t;
	return -1;
}

Lucas

首先可以观察得出 ，所以显然有

设

lucas定理也可以表示为一下形式

中 的系数显然为

可得 的系数为

因为 若 则 无法被表示系数为

得证。

1
int lucas(int n,int m,int mod){if(n==0&&m==0)return 1;return lucas(n/mod,m/mod,mod)*1LL*C(n%mod,m%mod,mod)%mod;}

Exlucas

这里的 可能不是质数

根据唯一分解定理

若 可以直接对于每个质数使用lucas求解，再使用CRT合并。

否则需要求解子问题 再使用CRT合并。

只需要解决 ，其中 是 包含的 的次数。

考虑将 分组，

可以发现后面是模 是循环的，可以通过预处理求出

前面的 是一个子问题递归求解即可。

复杂度 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
ll mul(ll a,ll b,ll mod)
{
	ll tmp=a*b-(ll)((long double)a/mod*b+0.5)*mod;
	return tmp<0?tmp+mod:tmp;
}
ll Pow(ll a,ll k,ll mod)
{
	ll ret=1;
	while(k){if(k&1)ret=ret*1LL*a%mod;k>>=1,a=a*1LL*a%mod;}
	return ret;
}
ll fact(ll n,ll p,ll pk)
{
	if(!n)return 1;
    ll ans=1;
    for(ll i=1;i<pk;i++)if(i%p)ans=ans*1LL*i%pk;
    ans=Pow(ans,n/pk,pk);
    for(ll i=1;i<=n%pk;i++)if(i%p)ans=ans*1LL*i%pk;
    return ans*1LL*fact(n/p,p,pk)%pk;
}
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
void exgcd(ll a,ll b,ll &x,ll &y)
{
	if(b==0)return x=1,y=0,void();
	exgcd(b,a%b,x,y);
	ll t=x;x=y,y=t-y*(a/b);
}
ll inv(ll x,ll p)
{
	ll t1,t2;
	exgcd(x,p,t1,t2);
	return (t1%p+p)%p;
}
ll getpower(ll x,ll p)
{
	ll c=0;
	while(x)c+=x/p,x/=p;
	return c;
}
ll C(ll n,ll m,ll p,ll pk)
{
	ll pa=getpower(n,p)-getpower(m,p)-getpower(n-m,p);
	return Pow(p,pa,pk)*1LL*fact(n,p,pk)%pk*inv(fact(m,p,pk),pk)%pk*inv(fact(n-m,p,pk),pk)%pk;
}
ll ExCRT(ll *a,ll *p,const ll&n)
{
	ll A=a[1],M=p[1],k1,k2;
	for(ll i=2;i<=n;i++)
	{
		ll c=((a[i]-A)%p[i]+p[i])%p[i],g=gcd(M,p[i]);
		if(c%g)return -1;
		exgcd(M,p[i],k1,k2);
		k1=mul(k1,c/g,p[i]);
		ll now=M/g*p[i];
		A=(mul(k1,M,now)+A)%now,M=now;
	}
	return (A%M+M)%M;
}
ll a[30],b[30];
ll cnt;
ll exlucas(ll n,ll m,ll p)
{
	cnt=0;
	ll lim=sqrt(p);
	for(ll i=2;i<lim;i++)
	{
		if(p%i)continue;
		ll t=1;while(p%i==0)p/=i,t*=i;
		a[++cnt]=t,b[cnt]=C(n,m,i,t);
	}
	if(p>1)a[++cnt]=p,b[cnt]=C(n,m,p,p);
	return ExCRT(b,a,cnt);
}

二次剩余

无特殊说明 均为奇素数。

对于 ，若存在 ，满足

则称 为模 意义下的二次剩余

勒让德符号

在模意义下是二次剩余在模意义下是非二次剩余

欧拉判别准则

证明：

若 为模 意义下的二次剩余，即存在 ，由费马小定理可得 ，显然 。

若 为模 意义下的非二次剩余，假设存在 ，此时 ，由费马小定理可知 不存在

若 ，显然成立。

Cipolla

定理1

定理2

的二次剩余与非二次剩余的个数均为 （不考虑 的情况下）。

证明：

我们只考虑所有的 ，假设有 且 ，则 ，即 ，若，则 ，故 ，也就是不同的 共有 个，即二次剩余有个 个。

算法流程
  1. 判断给定的数 是否是二次剩余。
  2. 随机一个 ，使其满足 是二次非剩余。
  3. 令 ，取 为其中一个可行解。

证明：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
int mod,w;
struct Complex
{
	int a,b;
	Complex(int A=0,int B=0){a=A,b=B;}
	Complex operator * (const Complex &t)const
	{
		return Complex((a*1LL*t.a+b*1LL*t.b%mod*w)%mod,(a*1LL*t.b+b*1LL*t.a)%mod);
	}
}; 
int Pow(int a,long long k)
{
	int ret=1;
	while(k){if(k&1)ret=ret*1LL*a%mod;k>>=1,a=a*1LL*a%mod;}
	return ret;
}
int Pow(Complex a,long long k)
{
	Complex ret=Complex(1,0);
	while(k){if(k&1)ret=ret*a;k>>=1,a=a*a;}
	return ret.a;
}
int Sqrt(int x)
{
	if(!x)return 0;
	if(Pow(x,(mod-1)>>1)==mod-1)return -1;
	while(true)
	{
		int a=rand()*1LL*rand()%mod;
		w=(a*1LL*a%mod+mod-x)%mod;
		if(Pow(w,(mod-1)>>1)==mod-1)return Pow(Complex(a,1),(mod+1)>>1);
	}
}

素数分解相关

Miller-Rabin

判断 是否为素数（）

根据费马小定理，当 且 为质数时，

而当 不是质数时，不一定成立

因此就有了一个想法：随机一些 ，对每一个 验证 ，如果有一个错误那么它一定是合数

但是，存在这样一类合数 ，对于 上面的式子都成立，比如

二次探测定理

如果 ，那么 或者 ，这里 是质数。

可以发现对于合数，有一定概率不成立。

在检验时，如果 且 ，则下一步计算 ,如果算出来不是 或 ，那么 是合数，如果算出来是 并且 ，则继续计算 ，直到出现奇数或者算出 。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
ll mul(ll a,ll b,ll mod)
{
	ll tmp=a*b-(ll)((long double)a/mod*b+0.5)*mod;
	return tmp<0?tmp+mod:tmp;
}
ll Pow(ll a,ll k,ll mod)
{
	ll ret=1;
	while(k){if(k&1)ret=mul(ret,a,mod);k>>=1,a=mul(a,a,mod);}
	return ret;
}
const int tprim[10]={0,2,3,5,7,11,13,17,19,23};
bool Miller_Rabin(ll x)
{
    if(x<2)return 0;
    ll k=x-1,cnt=0;
    while(!(k&1))k>>=1,++cnt;
    for(int i=1,j;i<=9;i++)
	{
		if(tprim[i]==x)return 1;
        if(Pow(tprim[i],x-1,x)!=1)return 0; 
        ll now=Pow(tprim[i],k,x);
        if(now==1||now==x-1)continue;
        now=mul(now,now,x);
        for(j=1;j<=cnt;++j,now=mul(now,now,x))if(now==x-1)break;	
        if(j>cnt)return 0;
    }
    return 1;
}

递归改成递推，这样复杂度是

莫比乌斯反演

莫比乌斯反演函数

其余情况

性质

莫比乌斯反演定理

若函数 满足

则存在

证明：

考虑 的系数可以很简单地得到证明。

推论：

若函数 满足

则存在

证明：

考虑每个 的系数被贡献时满足 ，显然 是 的约数，被贡献的系数值为 。

这个推论可以很简单地用于求解类似 问题

一些例题

基础的套路不在过多阐述。

Problem 1.

设 满足 （满足 的 的和）

显然

由莫比乌斯反演的性质可得

枚举

这里可以直接数论分块 解决。

继续将这个式子变形， 考虑枚举

可以发现后面的 是一个积性函数，可以线性筛出来。

这样对于单次询问复杂度

Problem 2.

为 的约数个数

有一个结论

显然每个质因子 的贡献是独立的， 为每种质因子贡献的乘积，设 可以发现左边贡献为 ，右边贡献为 （），所以得证。


按照套路考虑枚举 ，显然此时 对 有贡献，仅当 是 的倍数。

所以 为 的倍数， 为 的倍数。

考虑设 为

显然 不难发现

考虑线性筛筛出

这个显然可以数论分块，总复杂度

总结

一般来说莫比乌斯反演的题目，主要考虑化为 ，通过交换枚举顺序构造积性函数，或者使用数论分块求解，最重要的是抓住运算的性质和熟练掌握反演公式。

Min_25 筛

规定 是所有质数组成的集合，若无特殊说明 为质数。

其中 为积性函数，并且 能够快速计算， 能够表示为一个简单多项式。

筛质数的答案

首先需要对每个 求出

首先线性筛出 以内的所有质数， 表示第 个质数。

因为 能够写成一个简单多项式

所以本质上来说需要求解的是

考虑构造一个函数 满足

表示 的最小质因子。

其实不难发现 相当于运行了 次埃氏筛后，没有被筛掉的数和质数的 的和。

考虑如何求 ，发现这个转移可以分类讨论

考虑第 次筛质数，显然筛掉的是最小质因子大于 的数，而满足这个条件数必然大于

所以如果 ，这次筛是不会筛掉任何数的，显然 。

否则 ，显然会筛掉 的所有倍数（最小质因子为 ）。

不难发现 是一个完全积性函数，可以将需要删掉的 的倍数表示为 (，且 的最小质因子大于等于 )，不难发现 为所有合法的 的和加上 ，而 恰好等于 ，所以 就是所有合法的 的和，而 是完全积性函数，只需要乘上 就是需要删掉的数。

边界条件 显然 较小的情况可以直接手推公式， 较大可以考虑插值（ 大了时间复杂度爆炸）。

筛所有数的答案

如果分质数与合数讨论显然可以得到

前面一部分即为所有质数的贡献，而后面求合数的贡献，非常好理解就是暴力枚举 及其次幂然后求值，要注意的是 并不包含出 的函数值，所以还需要加上 。

总复杂度

关于实现方面的细节

有一个性质 所以我们只需要预处理所有的 就可以了。

在这里如果使用 map 储存复杂度显然会变大，考虑到 与 中必然有一个数 ，所以可以直接用两个数组储存编号。

运算的时候不需要记忆化，因为调用的所有 中没有相同的整数对 。

举个简单的例子

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
#include<bits/stdc++.h>
using namespace std;
namespace Min_25
{
	const int mod=1e9+7,inv2=(mod+1)/2,inv3=(mod+1)/3;
	const int maxn=1e6+10;
	int md(int x){return x>=mod?x-mod:x;}
	bool flag[maxn];
	int g1[maxn],g2[maxn];
	long long w[maxn];
	int id1[maxn],id2[maxn],m;
	int sum1[maxn],sum2[maxn],lim;	
	int prime[maxn],cnt=0;
	void make_prime(int n)
	{
		for(int i=2;i<=n;i++)
		{
			if(!flag[i])
			{
				prime[++cnt]=i;
				sum1[cnt]=(sum1[cnt-1]+i);
				sum2[cnt]=(sum2[cnt-1]+i*1LL*i)%mod;
			}
			for(int j=1;j<=cnt&&prime[j]*i<=n;j++)
			{
				flag[prime[j]*i]=1;
				if(i%prime[j]==0)break;
			}
		}
	}
	long long n;
	void G()
	{
		m=0;
		for(long long i=1,j;i<=n;i=j+1)
		{
			w[++m]=n/i,j=n/w[m];
			if(w[m]<=lim)id1[w[m]]=m;
			else id2[n/w[m]]=m;
			int t=w[m]%mod;
			g1[m]=(t+1)*1LL*t/2%mod-1;
			g2[m]=(t+1)*1LL*t/2%mod*(t*2+1)%mod*inv3%mod-1;
		}
		for(int i=1;i<=cnt;i++)
		for(int j=1;j<=m&&prime[i]*1LL*prime[i]<=w[j];j++)
		{
			int lst=w[j]/prime[i]<=lim?id1[w[j]/prime[i]]:id2[n/(w[j]/prime[i])];
			g1[j]=md(g1[j]+mod-(g1[lst]-sum1[i-1]+mod)*1LL*prime[i]%mod);
			g2[j]=md(g2[j]+mod-(g2[lst]-sum2[i-1]+mod)*1LL*prime[i]%mod*prime[i]%mod);
		}
	}
	int S(long long a,int j)
	{
		if(prime[j]>=a)return 0;
		int id=a<=lim?id1[a]:id2[n/a];
		int ans=md(md(g2[id]+mod-g1[id])+mod-md(sum2[j]+mod-sum1[j]));
		for(int k=j+1;prime[k]*1LL*prime[k]<=a&&k<=cnt;k++)
		{
			long long pe=prime[k];
			for(int e=1;pe<=a;e++,pe*=prime[k])
			{
				long long tmp=pe%mod;
				ans=(ans+tmp*(tmp-1)%mod*(S(a/pe,k)+(e!=1)))%mod;
			}
		}
		return ans;
	}	
	int Sum(long long N)
	{
		n=N;
		lim=sqrt(n);
		make_prime(lim);
		G();
		return (S(n,0)+1)%mod;
	}
}
long long n;
int main()
{
	scanf("%lld",&n);
	printf("%d\n",Min_25::Sum(n));
}

杜教筛

常见积性函数

  • 莫比乌斯反演函数。
  • 欧拉函数 。
  • 约数个数函数
  • 约数个数和函数
  • 恒等函数函数值恒为
  • 元函数
  • 单位函数

狄利克雷卷积

满足交换律，结合律，分配律。

不难发现存在 。

杜教筛

需要计算积性函数 的前缀和。

记

考虑寻找到两个积性函数使得 ，那么有

考虑将 对应的项提出来

也就是说只要 前缀和很好求，那么这个式子就可以使用整除分块完成计算。

例子

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
int smu[maxn];
long long sphi[maxn];
int mu[maxn],phi[maxn];
int prime[maxn],cnt=0,flag[maxn];
void make_prime(int n)
{
	mu[1]=1,phi[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(!flag[i])
		{
			prime[++cnt]=i;
			mu[i]=-1,phi[i]=i-1;
		}
		for(int j=1;j<=cnt&&prime[j]*i<=n;j++)
		{
			flag[prime[j]*i]=1;
			if(i%prime[j]==0)
			{
				mu[prime[j]*i]=0;
				phi[prime[j]*i]=phi[i]*prime[j];
				break;
			}
			mu[prime[j]*i]=-mu[i];
			phi[prime[j]*i]=phi[i]*phi[prime[j]];
		}
	}
	for(int i=1;i<=n;i++)smu[i]=smu[i-1]+mu[i],sphi[i]=sphi[i-1]+phi[i];
}
unordered_map<long long,long long>phimp;
unordered_map<long long,int>mump;
long long Sphi(long long n)
{
	if(n<=5e6)return sphi[n];
	if(phimp.count(n))return phimp[n];
	long long ans=n*(n+1)/2;
	for(long long i=2,j;i<=n;i=j+1)
	{
		j=n/(n/i);
		ans-=(j-i+1)*Sphi(n/i);
	}
	return phimp[n]=ans;
}
int Smu(long long n)
{
	if(n<=5e6)return smu[n];
	if(mump.count(n))return mump[n];
	int ans=1;
	for(long long i=2,j;i<=n;i=j+1)
	{
		j=n/(n/i);
		ans-=(j-i+1)*Smu(n/i);
	}
	return mump[n]=ans;
}
  • Title: 数论算法
  • Author: Reaepita
  • Created at: 2022-09-09 19:47:59
  • Updated at: 2023-08-06 15:21:37
  • Link: https://harrybh.github.io/2022/09/09/数论总结/
  • License: This work is licensed under CC BY-NC-SA 4.0.
 Comments
On this page
数论算法
  1. 数论算法
    1. 快速乘
    2. 同余系相关
      1. Exgcd
      2. ExCRT
      3. ExBSGS
      4. Lucas
      5. Exlucas
      6. 二次剩余
    3. 素数分解相关
      1. Miller-Rabin
    4. 莫比乌斯反演
      1. 莫比乌斯反演函数
      2. 莫比乌斯反演定理
      3. 一些例题
      4. 总结
    5. Min_25 筛
      1. 筛质数的答案
      2. 筛所有数的答案
      3. 关于实现方面的细节
    6. 杜教筛
      1. 常见积性函数
      2. 狄利克雷卷积
      3. 杜教筛
      4. 例子