数论难点选讲
Miller_Rabin素性测试
我们首先看这样一个很简单的问题：判定正整数\(n\)是否为素数
最简单的做法就是枚举\(2\)到\(n\)的所有数，看是否有数是\(n\)的因数，时间复杂度\(O(n)\)
稍微优化一下发现只要枚举\(2\)到\(\sqrt{n}\)中的数就可以了
然后发现数据范围\(n\leq 10^{18}\)，时间复杂度直接就死掉了QAQ
我们就要考虑新的方法了
首先引入两个定理
1、费马小定理
如果\(p\)是素数，且\(gcd(a,p)=1\)，那么\(a^{p-1}\equiv 1(mod \ n)\)
证明什么的你随便找本数论书自己翻一下
注意它的逆定理不一定成立（或者说是它的逆定理在大多数情况下都成立）
2、二次探测定理（其实这也没有一个准确的名字）
如果\(p\)是奇素数，\(x<p\)，且\(x^2\equiv1(mod\ p)\)，那么\(x=1\)或\(x=p-1\)
证明：由同余式知\(x^2-1\equiv0(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}\text%n\)，如果这个值不为1那么\(n\)一定是合数（利用费马小定理）
3、遍历\(i\)，使得\(1\leq i \leq q\)，如果\(2^i*m\text%n=1\)并且\(a^{i-1}*m\text%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\n");return;}
//cout <<1;
int l=n-1,r=n+1;
while (!miller_rabin(l)) l--;
while (!miller_rabin(r)) r++;
printf("%d\n",r-l);
}
signed main()
{
n=read();
while (n)
{
work();
n=read();
}
return 0;
}
BSGS相关
BSGS
给定\(a,b,p\)，求最小的\(x\)使得\(a^x\equiv b(mod\ p)\)，保证\(p\)是质数
考虑这样的一个数\(m\)，使得\(x=qm-r(q\in[1...p/m],r\in[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 \rceil\)时间复杂度最优，为\(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\nmid b\)，那么就只有在\(y=1,x=0\)的情况下有解，其他情况均无解
就这样一直除下去，直到\(d=1\)，时停止。这时你得到了一个\(d\)的序列\(d_1,d_2,\cdots,d_t\)
以及一个变形之后的同余式：\(a^{x-t}\frac{x^t}{\prod d_i}\equiv\frac{b}{\prod d_i}(mod\ \frac{p}{\prod d_i})\)
好像直接BSGS就可以了？
是的，但是要先特判掉\(x\in[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\n");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\n");ok=0;break;}
b/=d;p/=d;k=(k*a/d)%p;t++;
if (k==b) {printf("%d\n",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\n",m*i-mp[sum]+t);ok=0;break;}
}
if (ok) printf("No Solution\n");
}
}
return 0;
}
原根
记住：\(1004535809\)和\(998244353\)的原根为\(3\)，\(1000000007\)的原根为\(5\)
相关定义
阶：若\(gcd(a,p)=1\)，则最小的正整数\(n\)使得\(a^n\equiv 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 i\leq m-2)\)取遍\(1\)到\(p-1\)
求法
暴力的想法是\(g\)从\(2\)开始枚举，检验当前\(g\)模\(p\)的阶是否为\(\varphi(p)\)，检验阶的方法就是让\(i\)从\(1\)取到\(p-1\)，判断是否出现了\(g^i\equiv 1(mod\ p)\)
考虑优化，前面对\(g\)的枚举似乎不好优化，考虑检验
结论：若对\(\forall i|\varphi(p)\)，都没有\(g^i\equiv1(mod\ p)\)，那么\(\forall i,0\leq i\leq \varphi(p)\)，均无\(g^i\equiv 1(mod\ p)\)（\(i\neq \varphi(p)\)）
证明：首先有\(g^{\varphi(p)}\equiv 1(mod\ p)\)
考虑反证，假设存在一个最小的数\(c\)，满足\(g^c\equiv 1(mod\ p)\)，那么有条件值\(c<\varphi(p)\)且\(c\)不是\(\varphi(p)\)的因数）
记\(d=\varphi(p)-c\)，那么有\(g^d\equiv 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（d\in 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\)
其中满足\(gcd(b_i,b_j)=1,\forall i,j ,i\neq j\)
方法：
1）记\(B=b_1b_2\cdots b_m\)
2）求解\(m\)个同余方程组\(\frac{B}{b_i}t_i\equiv1(mod\ b_i)\)
3）\(ans=\sum_{i=1}^m\frac{B}{b_i}t_ia_i\)，同时\(ans-=B*k\)使得\(k\)最大且\(ans\)不为负
证明（假的）：
对于第\(i\)个同余方程组，因为\(b_j|B\)且\(gcd(b_i,b_j)=1\)，所以这一组同余方程得到的解不会影响到其它方程的解，
又因为最后减去的数是所有的\(b\)的倍数，所以也不会影响到同余式的性质
#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;
}
扩展型
非常抱歉，这篇文章鸽了
