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)
其中满足(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)的倍数，所以也不会影响到同余式的性质
#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;
}
扩展型
非常抱歉，这篇文章鸽了