数论难点选讲

Miller_Rabin素性测试

我们首先看这样一个很简单的问题：判定正整数(n)是否为素数

最简单的做法就是枚举(2)到(n)的所有数，看是否有数是(n)的因数，时间复杂度(O(n))

稍微优化一下发现只要枚举(2)到(sqrt{n})中的数就可以了

然后发现数据范围(nleq 10^{18})，时间复杂度直接就死掉了QAQ

我们就要考虑新的方法了

首先引入两个定理

1、费马小定理

如果(p)是素数，且(gcd(a,p)=1)，那么(a^{p-1}equiv 1(mod n))

证明什么的你随便找本数论书自己翻一下

注意它的逆定理不一定成立（或者说是它的逆定理在大多数情况下都成立）

2、二次探测定理（其实这也没有一个准确的名字）

如果(p)是奇素数，(x<p)，且(x^2equiv1(mod p))，那么(x=1)或(x=p-1)

证明：由同余式知(x^2-1equiv0(mod p))，即(p|(x+1)(x-1))

​ 又由(p)是素数知(p|x-1)或(p|x+1)，解得(x=1)或(x=p-1)

诶等等zzr没事给证明干嘛？zzr不是最讨厌证明了吗

由上面很简单的证明过程我们可以发现，(x=1)和(x=p-1)这两个解其实是对所有的(p)都成立的

即无论(p)取什么值(x)取上面两个值是一定可以的

但是当(p)是一个合数的时候，此时原同余方程的解(x)就不只上面这两个了，而是会有多个

换一句话说：如果上面的(x)取到了1和(p-1)以外的数，就说明(p)不是一个素数了

我们主要利用上面两个性质来进行素数判定

1、取(2^q*m=n-1)（(q,m)均为正整数且(m)为奇数），同时任意取小于(n)的正整数(a)

2、求出(a^{n-1} ext%n)，如果这个值不为1那么(n)一定是合数（利用费马小定理）

3、遍历(i)，使得(1leq i leq q)，如果(2^i*m ext%n=1)并且(a^{i-1}*m ext%n!=1或n-1)，那么由二次探测定理就知道原同余方程出现一个特殊解，说明(n)不是一个素数

上面的方法有一个小问题：由于费马小定理的逆定理不一定成立（在大多数情况下成立），因此有时我们会对(n)进行误判，具体的，每做一次发生误判的概率是(frac{1}{4})

解决的方法在上面的解法中也有体现：换用不同的(a)，多进行几次即可

好了上面就是完整的miller-rabin测试了

一道例题：poj3518Prime Gap

题意：两个相邻的素数的差值叫做Prime Gap。输入一个K，求K两端的素数之差，如果K本身是一个素数，输出0；

分析：其实数据很小你直接筛一下也可以

​ 或者你直接暴力寻找当前这个数相邻的数是否是质数，两端分别记录第一次找到的质数即可

#include<iostream>
#include<string.h>
#include<string>
#include<stdio.h>
#include<stdlib.h> 
#include<algorithm>
#include<vector>
#include<queue>
#include<map>
using namespace std;
#define int long long
int n;

int read()
{
	int x=0,f=1;char ch=getchar();
	while ((ch<'0') || (ch>'9')) {if (ch=='-') f=-1;ch=getchar();}
	while ((ch>='0') && (ch<='9')) {x=x*10+(ch-'0');ch=getchar();}
	return x*f;
}

int mul(int x,int y,int n)
{	
	x%=n;y%=n;
	int ans=0,sum=x;
	while (y)
	{
		int tmp=y%2;y/=2;
		if (tmp) ans=(ans+sum)%n;
		sum=(sum+sum)%n;
	}
	return ans;
}

int qpow(int x,int y,int n)
{
	int ans=1,sum=x;
	while (y)
	{
		int tmp=y%2;y/=2;
		if (tmp) ans=mul(ans,sum,n);
		sum=mul(sum,sum,n);
	}
	return ans;
}

bool prime(int m,int q,int a,int n)
{
	int now=qpow(a,m,n);
	if ((now==1) || (now==n-1)) return 1;
	int i;
	for (i=1;i<=q;i++)
	{
		int x=mul(now,now,n);
		if ((x==1) && (now!=1) && (now!=n-1)) return 0;
		now=x;
	}
	if (now!=1) return 0;//其实这里是将费马小定理的检测放在了最后，省去再做一次快速幂
	return 1;
}

bool miller_rabin(int x)
{
	if (x==2) return 1;
	if ((x<2) || (x%2==0)) return 0;
	int num=x-1,tim=0;
	while ((num) && (num%2==0)) {num/=2;tim++;}
	//cout << num << " " <<tim << endl;
	int i;
	for (i=1;i<=10;i++)//一般都会进行20次左右，不过数据范围小对吧2333
	{
		int a=rand()%(x-1)+1;
		if (!prime(num,tim,a,x)) return 0;
	}
	return 1;
}

void work()
{
	if (miller_rabin(n)) {printf("0
");return;}
	//cout <<1;
	int l=n-1,r=n+1;
	while (!miller_rabin(l)) l--;
	while (!miller_rabin(r)) r++;
	printf("%d
",r-l);
}

signed main()
{
	n=read();
	while (n)
    {
		work();
		n=read();
	}
	return 0;
}

BSGS相关

BSGS

给定(a,b,p)，求最小的(x)使得(a^xequiv b(mod p))，保证(p)是质数

考虑这样的一个数(m)，使得(x=qm-r(qin[1...p/m],rin[0...m)))

那么原式则为(a^{qm-r}equiv b(mod p))，即(a^{qm}equiv b·a^r(mod p))

我们暴力枚举(r)，将(b·a^r)的值存入一个map中

接着枚举(q)，计算出对应的(a^{qm})，如果这个值在map中存在的话就直接输出结果，否则继续搜索直到无解

时间复杂度(O(frac{p}{m}·m))，取(m=lceil sqrt p ceil)时间复杂度最优，为(O(sqrt(p)))

模板题：[TJOI2007]可爱的质数

#include<iostream>
#include<string.h>
#include<string>
#include<stdio.h>
#include<algorithm>
#include<math.h>
#include<vector>
#include<queue>
#include<map>
using namespace std;
const int maxd=1000000007,N=100000;
const double pi=acos(-1.0);
typedef long long ll;
ll a,b,c;
map<ll,int> mp;

int read()
{
	int x=0,f=1;char ch=getchar();
	while ((ch<'0') || (ch>'9')) {if (ch=='-') f=-1;ch=getchar();}
	while ((ch>='0') && (ch<='9')) {x=x*10+(ch-'0');ch=getchar();}
	return x*f;
}

ll qpow(ll x,ll y)
{
	ll ans=1,sum=x;
	while (y)
	{
	    int tmp=y%2;y/=2;
	    if (tmp) ans=(ans*sum)%c;
	    sum=(sum*sum)%c;
	}
	return ans;
}

int main()
{
	c=read();a=read();b=read();
	int i,m=sqrt(c)+1;ll sum=b;
	for (i=0;i<m;i++)
	{
		mp[sum]=i;
		sum=(sum*a)%c;
	}
	ll tmp=qpow(a,m);sum=1;
	for (i=1;i<=m;i++)
	{
		sum=(sum*tmp)%c;
		if (mp.count(sum)) {printf("%d",i*m-mp[sum]);return 0;}
	}
	printf("no solution");
	return 0;
}

EXBSGS

还是上面那个问题，如果此时(p)不是质数呢？

上面的分块做法似乎就不是很可行了因为保证了(gcd(a,p)=1)所以对于不同的(x)，(a^x)模(p)的值两两不同

对于这个问题我们考虑转化为BSGS，我们在同余式两边同时除以(gcd(a,p))，记作(d)

注意此时如果(d mid b)，那么就只有在(y=1,x=0)的情况下有解，其他情况均无解

就这样一直除下去，直到(d=1)，时停止。这时你得到了一个(d)的序列(d_1,d_2,cdots,d_t)

以及一个变形之后的同余式：(a^{x-t}frac{x^t}{prod d_i}equivfrac{b}{prod d_i}(mod frac{p}{prod d_i}))

好像直接BSGS就可以了？

是的，但是要先特判掉(xin[0..t])的情况

模板题：SPOJ MOD - Power Modulo Inverted

#include<iostream>
#include<string.h>
#include<string>
#include<stdio.h>
#include<algorithm>
#include<math.h>
#include<vector>
#include<queue>
#include<map>
using namespace std;
const int maxd=1000000007,N=100000;
const double pi=acos(-1.0);
typedef long long ll;
map<ll,int> mp;
ll a,b,p;

int read()
{
	int x=0,f=1;char ch=getchar();
	while ((ch<'0') || (ch>'9')) {if (ch=='-') f=-1;ch=getchar();}
	while ((ch>='0') && (ch<='9')) {x=x*10+(ch-'0');ch=getchar();}
	return x*f;
}

ll qpow(ll x,ll y)
{
	ll ans=1,sum=x;
	while (y)
	{
		int tmp=y%2;y/=2;
		if (tmp) ans=(ans*sum)%p;
		sum=(sum*sum)%p;
	}
	return ans;
}

ll gcd(ll x,ll y)
{
	if (!y) return x;else return gcd(y,x%y);
}

int main()
{
	while (1)
	{
		a=read();p=read();b=read();
		if ((!a) && (!p) && (!b)) break;
		a%=p;b%=p;
		if (b==1) {printf("0
");continue;}
		ll d=gcd(a,p),k=1,t=0;bool ok=1;
		//cout << "d=" << d << endl;
		while (d!=1)
		{
		    if (b%d) {printf("No Solution
");ok=0;break;} 
		    b/=d;p/=d;k=(k*a/d)%p;t++;
		    if (k==b) {printf("%d
",t);ok=0;break;}
		    d=gcd(a,p);
		}
		if (ok)
		{
			mp.clear();int i;ll sum=b;
			ll m=sqrt(p)+1;
			for (i=0;i<m;i++) 
			{
				mp[sum]=i;
				//cout << sum << " ";
				sum=(sum*a)%p;
			}
			//cout << endl;
			//cout << t << endl;
			ll tmp=qpow(a,m);sum=k;
			for (i=1;i<=m;i++)
			{
				sum=(sum*tmp)%p;
				if (mp.count(sum)) {printf("%lld
",m*i-mp[sum]+t);ok=0;break;}
			}
			if (ok) printf("No Solution
");
	    }
	}
	return 0;
}

原根

记住：(1004535809)和(998244353)的原根为(3)，(1000000007)的原根为(5)

相关定义

阶：若(gcd(a,p)=1)，则最小的正整数(n)使得(a^nequiv 1)被称作(a)模(p)的阶

因为(a^{varphi(p)}equiv 1)，所以一定有(a)模(p)的阶(n|varphi(p))

若(a)模(p)的阶为(varphi(p))，那么(a)为模(p)的一个原根

性质

1、只有当(a=2,4,p^a,2*p^a)时，(a)具有原根

2、记模(m)的原根为(g)，则(g^i(0 leq ileq m-2))取遍(1)到(p-1)

求法

暴力的想法是(g)从(2)开始枚举，检验当前(g)模(p)的阶是否为(varphi(p))，检验阶的方法就是让(i)从(1)取到(p-1)，判断是否出现了(g^iequiv 1(mod p))

考虑优化，前面对(g)的枚举似乎不好优化，考虑检验

结论：若对(forall i|varphi(p))，都没有(g^iequiv1(mod p))，那么(forall i,0leq ileq varphi(p))，均无(g^iequiv 1(mod p))（(i eq varphi(p))）

证明：首先有(g^{varphi(p)}equiv 1(mod p))

考虑反证，假设存在一个最小的数(c)，满足(g^cequiv 1(mod p))，那么有条件值(c<varphi(p))且(c)不是(varphi(p))的因数）

记(d=varphi(p)-c)，那么有(g^dequiv 1(mod p))，由假设知(d>c)

那么我们又可以推出(g^{d-c}equiv 1(mod p))

由更相减损术知(g^{gcd(c,d)}equiv 1(mod p))，因为(gcd(c,d)leq c)，又由假设知(gcd(c,d)=c)，也就是(c|d)

设(d=kc（din N_+))，则(varphi(p)=c+d=(k+1)c)，推得(c|varphi(p))，与条件矛盾，QED

于是我们在验证(g)的时候只需要枚举(varphi(p))的因数进行判断即可

时间已经很优秀了，但是我们还可以继续优化

将(varphi(p))分解质因数(varphi(p)=p_1^{alpha_1}p_2^{alpha_2}cdots p_m^{alpha_m})，于是只需要检验(g^{frac{varphi(p)}{p_i}}mod p)即可

因为此时我们用来检验的幂指数都是(varphi(p))中小于(p)的因数，它至少比(p)少了一个质因子，并且(1)的任何正整数次幂都是(1)

于是就有了下面一份代码（模板：51nod1135）

#include<iostream>
#include<string.h>
#include<string>
#include<stdio.h>
#include<algorithm>
#include<math.h>
#include<vector>
#include<queue>
#include<map>
#include<set>
using namespace std;
#define lowbit(x) (x)&(-x)
#define rep(i,a,b) for (int i=a;i<=b;i++)
#define per(i,a,b) for (int i=a;i>=b;i--)
#define maxd 1000000007
typedef long long ll;
const int N=100000;
const double pi=acos(-1.0);
ll m,cnt=0,q[1001000];

int read()
{
	int x=0,f=1;char ch=getchar();
	while ((ch<'0') || (ch>'9')) {if (ch=='-') f=-1;ch=getchar();}
	while ((ch>='0') && (ch<='9')) {x=x*10+(ch-'0');ch=getchar();}
	return x*f;
}

ll qpow(ll x,ll y,ll p)
{
	ll ans=1;
	while (y)
	{
		if (y&1) ans=(ans*x)%p;
		x=(x*x)%p;
		y>>=1;
	}
	return ans;
}

ll get_phi(ll x)
{
	ll ans=x,i;
	for (i=2;i*i<=x;i++)
	{
		if (x%i) continue;
		ans=ans/i*(i-1);
		while (x%i==0) x/=i;
	}
	if (x>1) ans=ans/x*(x-1);
	return ans;
}

ll get_g(ll m)
{
	ll i,phi=get_phi(m);
	cnt=0;ll rest=phi;
	for (i=2;i*i<=m;i++)
	{
		if (phi%i) continue;
		q[++cnt]=phi/i;
		while (rest%i==0) rest/=i;
	}
	if (rest>1) q[++cnt]=phi/rest;
	ll g=2;
	while (1)
	{
		if (qpow(g,phi,m)!=1) {g++;continue;}
		bool flag=1;
		rep(i,1,cnt)
			if (qpow(g,q[i],m)==1) {g++;flag=0;break;}
		if (flag) return g;
	}
}

int main()
{
	m=read();
	printf("%lld",get_g(m));
	return 0;
}

作用

化乘法为加法，从而改变转移方程的形式，便于使用基于加法转移的优化操作

比如矩阵快速幂（CF1106F），FFT（SDOI序列统计）

中国剩余定理（CRT）

普通型

求解符合下列同余方程式的最小的(n)

[egin{cases} nequiv a_1(mod b_1)\ nequiv a_2(mod b_2)\ cdots\ nequiv a_m(mod b_m) end{cases} ]

其中满足(gcd(b_i,b_j)=1,forall i,j ,i eq j)

方法：

1）记(B=b_1b_2cdots b_m)

2）求解(m)个同余方程组(frac{B}{b_i}t_iequiv1(mod b_i))

3）(ans=sum_{i=1}^mfrac{B}{b_i}t_ia_i)，同时(ans-=B*k)使得(k)最大且(ans)不为负

证明（假的）：

对于第(i)个同余方程组，因为(b_j|B)且(gcd(b_i,b_j)=1)，所以这一组同余方程得到的解不会影响到其它方程的解，

又因为最后减去的数是所有的(b)的倍数，所以也不会影响到同余式的性质

代码luogu 3868

#include<iostream>
#include<string.h>
#include<string>
#include<stdio.h>
#include<algorithm>
#include<math.h>
#include<vector>
#include<queue>
#include<map>
#include<set>
using namespace std;
#define lowbit(x) (x)&(-x)
#define sqr(x) (x)*(x)
#define fir first
#define sec second
#define rep(i,a,b) for (register int i=a;i<=b;i++)
#define per(i,a,b) for (register int i=a;i>=b;i--)
#define maxd 1000000007
#define eps 1e-6
typedef long long ll;
const int N=100000;
const double pi=acos(-1.0);
int n;
ll a[20],b[20],x[20],lcm;

ll read()
{
	ll x=0,f=1;char ch=getchar();
	while ((ch<'0') || (ch>'9')) {if (ch=='-') f=-1;ch=getchar();}
	while ((ch>='0') && (ch<='9')) {x=x*10+(ch-'0');ch=getchar();}
	return x*f;
}

void exgcd(ll a,ll b,ll &x,ll &y)
{
	if (b==0) {x=1;y=0;return;}
	else
	{
		exgcd(b,a%b,x,y);
		ll tmpx=x,tmpy=y;
		x=tmpy;y=tmpx-a/b*tmpy;
	}
}

ll mul(ll x,ll y)
{
	ll ans=0;
	while (y)
	{
		if (y&1) ans=(ans+x)%lcm;
		x=(x+x)%lcm;
		y>>=1;
	}
	return ans;
}

int main()
{
	n=read();
	rep(i,1,n) a[i]=read();
	rep(i,1,n) b[i]=read();
	rep(i,1,n) a[i]=(a[i]%b[i]+b[i])%b[i];
	lcm=1;ll ans=0,y;
	rep(i,1,n) lcm*=b[i];
	rep(i,1,n)
	{
		exgcd(lcm/b[i],b[i],x[i],y);
		x[i]=(x[i]%b[i]+b[i])%b[i];
		ans=(ans+mul(mul(lcm/b[i],x[i]),a[i]))%lcm;
	}
	ans=(ans+lcm)%lcm;
	printf("%lld",ans);
	return 0;
}

扩展型

非常抱歉，这篇文章鸽了