先记录一下一些概念和定理
同余:给定整数a,b,c,若用c不停的去除a和b最终所得余数一样,则称a和b对模c同余,记做a≡b (mod c),同余满足自反性,对称性,传递性
定理1:
若a≡b (mod c),对某个整数k有
a+k≡b+k (mod c)
a-k≡b-k (mod c)
ak≡bk (mod c)
定理2:
若a≡b (mod c),d≡e (mod c),有
ax+dy≡bx+ey (mod c) ,x,y为任意整数,即同余式可以相加
ad≡be (mod c) ,即同余式可以相乘
an≡bn (mod c),n>0
f(a)≡f(b) (mod c) f(x)为任一整系数多项式
定理3:
若a≡b (mod c),且d|c 则a≡b (mod d)
若a≡b (mod c),则gcd(a,c)=gcd(b,c)
a≡b (mod ci)(1<=i<=n) 同时成立 当且仅当a若a≡b (mod lcm(c1,c2,c3...cn))
若ad≡bd (mod c)且gcd(d,c)=e 则a≡b (mod c/e)
然后就是一元线性同余方程了,也就是形如 axΞb (mod c) 这样的方程
那么上面的方程可以写作ax+cy=b的二元一次方程,那很明显这个就可以用扩展欧几里得来解决。
设g=gcd(a,c),若b%g!=0则方程无解,否则a0=a/g,b0=b/g,c0=c/g,此时gcd(a0,c0)=1,
a0x+c0y=b0就可以用扩展欧几里得解出x,然后通解就是x+kb0
至于扩展欧几里得就不多重复了,直接挂两个博客
Yo!Hadilo!的欧几里得算法与扩展欧几里得算法_C++
镜外之主 的ACM数论之旅4---扩展欧几里德算法(欧几里德(・∀・)?是谁?)
然后就是一元线性方程组了,也就是一组一元线性方程组
x≡b1(mod c1)
x≡b2(mod c2)
x≡b3(mod c3)
....
x≡bn(mod cn)
那么我们就由前三组来看,
x=k1*c1+b1,x=k2*c2+b2,x=k3*c3+b3 ,k1,k2,k3都为任意整数。
前两个式子联立得,k1*c1+b1=k2*c2+b2,整理可得k2*c2-k1*c1=b1-b2,
此时扩展欧几里得可以解出一个最小非负整数解k2",或者直接无解了。
那么x有一个解为x2=k2"*c2+b2,而k2的通解为k2=k2"+m*c1/gcd(c1,c2), m为任意整数。
再引入第三个式子,就有k3*c3+b3=k2*c2+b2。
代入k2 也就是k3*c3+b3=k2"*c2+m*c1*c2/gcd(c1,c2)+b2,再代入x2整理就得,
k3*c3-m*c1*c2/gcd(c1,c2)=x2-b3,这时解出k3",就有x3,和k3通解
那么以此类推,联立后面的方程,最终解出的xn就是最终结果了。
Strange Way to Express Integers POJ - 2891
1 #include<cstdio> 2 typedef long long ll; 3 const int N=1e5+11; 4 ll bb[N],cc[N]; 5 ll exgcd(ll a,ll b,ll &x,ll &y){ 6 if(!b){ 7 x=1; 8 y=0; 9 return a; 10 } 11 ll g=exgcd(b,a%b,y,x); 12 y-=a/b*x; 13 return g; 14 } 15 ll linemod(int n){ 16 ll a,b,c,ans=bb[0],cp=cc[0],x,y,g; 17 for(int i=1;i<n;i++){ 18 a=cc[i];c=cp;b=ans-bb[i]; 19 g=exgcd(a,c,x,y); 20 if(b%g) return -1; 21 a/=g;b/=g;c/=g; 22 x=(x*b%c+c)%c; 23 ans=x*cc[i]+bb[i]; 24 cp*=a; 25 ans=(ans%cp+cp)%cp; 26 } 27 return ans; 28 } 29 int main(){ 30 int n; 31 while(~scanf("%d",&n)){ 32 for(int i=0;i<n;i++) scanf("%lld%lld",&cc[i],&bb[i]); 33 printf("%lld ",linemod(n)); 34 } 35 return 0; 36 }
还有几题 HDU1021 HDU 1573 HDU3579都可以练一练手,HDU国庆进不去,嘤嘤嘤。
然后多元线性同余方程的话就简单把书上的抄一下。
定义:形如a1x2+a2x2+a3x3+...+anxn+b≡0(mod c) 的同余方程成为多元线性同余方程。
性质:多元线性同余方程。a1x2+a2x2+a3x3+...+anxn+b≡0(mod c) 有解的充分必要条件是gcd(a1,a2,a3,...an,c) | b,若有解,解的个数为cn-1*gcd(a1,a2,a3...,an,c)
解法:令d=gcd(a1,a2,a3,...an,c) ,d1=gcd(a1,a2,a3,...an-1,c),那么gcd(d1,an)=d,由方程可知anxn+b≡0 (mod d1)
这样那方程分成了a1x2+a2x2+a3x3+...+an-1xn-1+b1d1≡0(mod c) (b1d1=anxn+b) 和anxn+b≡0 (mod d1)两个方程,两个都满足时即可一步步解出x
这个没啥题,就当扩展知识,了解一下。
然后是高次同余方程,我的书里说,高次同余方程比较复杂,所以本书仅讨论ax≡b(mod c) (0<=x<C)的问题
然后解决的方法是拔山盖世BSGS算法,全称是Baby Step Giant Step,小步大步算法,因为整个算法取决于一个分块的思想。
先放几个博客,因为在构造方式上的不同,所以在实现方法上有着要不要通过扩展欧几里得算法来求逆元的差异,我比较懒,就没有。
首先我们先来对于c是质数来进行讨论,c是质数的话,那么由费马小定理就有aphi(c)≡1(mod c),那么ax≡b(mod c) 有解的话,x肯定是在0~phi(c) phi(c)是c的欧拉函数哦
那么我们可以设x=i*n-j,n为√c,0<=j<n那么就有ai*n-j≡b(mod c) ,也就是ai*n≡b*aj(mod c),(设x=i*n+j的话,就得用扩展欧几里得求逆元了)
上面就是分块的思想,每个an为一块,这样我们先枚举j,将b*aj预处理保存在哈希表中,然后在枚举i,看相应的ai*n在哈希表有没有存在,存在的话,就可以得到一组(i,j),
代入x=i*n-j,即可得出x,总的时间复杂度就√c
而如果c不是质数,也是a和c不互质的话,那么此时就是想办法使得gcd(a,c)=1,这个在上面的扩展BSGS里详细的证明了,这里不重复了。
Discrete Logging POJ - 2417
Clever YPOJ - 3243
Mod Tree HDU - 2815
上面是三个板子题,可以练练手。
1 #include<cstdio> 2 #include<cmath> 3 #include<algorithm> 4 using namespace std; 5 typedef long long ll; 6 const int inf=1e9+7; 7 struct Hashmap{ 8 static const int Hash=999917,N=46340; 9 int num,head[Hash],a[N+11],b[N+11],ne[N+11]; 10 int top,sta[N+11]; 11 void clear(){ 12 num=0; 13 while(top) head[sta[top--]]=0; 14 } 15 void add(int x,int y){ 16 a[++num]=y; 17 ne[num]=head[x]; 18 b[num]=inf; 19 head[x]=num; 20 } 21 bool count(int y){ 22 int x=y%Hash; 23 for(int j=head[x];j;j=ne[j]) 24 if(y==a[j]) return true; 25 return false; 26 } 27 int &operator [](int y){ 28 int x=y%Hash; 29 for(int j=head[x];j;j=ne[j]) 30 if(y==a[j]) return b[j]; 31 add(x,y); 32 sta[++top]=x; 33 return b[num]; 34 } 35 }mmp; 36 ll bsgs(ll a,ll b,ll c){ 37 if(c==1) return !b ? !a : -1; 38 if(a%c==0) return !b ? 1 : -1; 39 if(b==1) return 0; 40 int sqc=ceil(sqrt(1.0*c)); 41 ll aa=1,ac=1,ab; 42 mmp.clear(); 43 for(int i=0;i<sqc;i++){ 44 ab=aa*b; 45 if(ab>=c) ab%=c; 46 mmp[ab]=i; 47 aa*=a; 48 if(aa>=c) aa%=c; 49 } 50 for(int i=1;i<=sqc;i++){ 51 ac*=aa; 52 if(ac>=c) ac%=c; 53 if(mmp.count(ac)) return i*sqc-mmp[ac]; 54 } 55 return -1; 56 } 57 int main(){ 58 ll x,y,z,ans; 59 while(~scanf("%lld%lld%lld",&z,&x,&y)){ 60 ans=bsgs(x,y,z); 61 if(ans==-1) printf("no solution "); 62 else printf("%lld ",ans); 63 } 64 return 0; 65 }
1 #include<cstdio> 2 #include<cmath> 3 #include<map> 4 #include<algorithm> 5 using namespace std; 6 typedef long long ll; 7 const int inf=1e9+7; 8 struct Hashmap{ 9 static const int Hash=999917,N=1e5+11; 10 int num,head[Hash],a[N+11],b[N+11],ne[N+11]; 11 int top,sta[N+11]; 12 void clear(){ 13 num=0; 14 while(top) head[sta[top--]]=0; 15 } 16 void add(int x,int y){ 17 a[++num]=y; 18 ne[num]=head[x]; 19 b[num]=inf; 20 head[x]=num; 21 } 22 bool count(int y){ 23 int x=y%Hash; 24 for(int j=head[x];j;j=ne[j]) 25 if(y==a[j]) return true; 26 return false; 27 } 28 int &operator [](int y){ 29 int x=y%Hash; 30 for(int j=head[x];j;j=ne[j]) 31 if(y==a[j]) return b[j]; 32 add(x,y); 33 sta[++top]=x; 34 return b[num]; 35 } 36 }mmp; 37 //map<ll,int> mmp; 38 ll gcd(ll a,ll b){ 39 return !b ? a : gcd(b,a%b); 40 } 41 ll bsgs(ll a,ll b,ll c){ 42 if(b>=c) return -1; 43 if(c==1) return !b ? !a : -1; 44 if(a%c==0) return !b ? 1 : -1; 45 if(b==1) return 0; 46 int sqc,d=0; 47 ll aa=1,ac=1,ab,g; 48 while((g=gcd(a,c))!=1){ 49 if(b%g) return -1; 50 b/=g;c/=g; 51 ac*=a/g; 52 if(ac>=c) ac%=c; 53 d++; 54 if(ac==b) return d; 55 } 56 mmp.clear(); 57 sqc=ceil(sqrt(1.0*c)); 58 for(int i=0;i<sqc;i++){ 59 ab=aa*b; 60 if(ab>=c) ab%=c; 61 mmp[ab]=i; 62 aa*=a; 63 if(aa>=c) aa%=c; 64 } 65 for(int i=1;i<=sqc;i++){ 66 ac*=aa; 67 if(ac>=c) ac%=c; 68 if(mmp.count(ac)) return sqc*i-mmp[ac]+d; 69 } 70 return -1; 71 } 72 int main(){ 73 ll x,y,z,ans; 74 while(~scanf("%lld%lld%lld",&x,&z,&y)){ 75 ans=bsgs(x,y,z); 76 if(ans==-1) printf("Orz,I can’t find D! "); 77 else printf("%lld ",ans); 78 } 79 return 0; 80 }
最后就是中国剩余定理了。直接上两博客吧。
镜外之主的ACM数论之旅9---中国剩余定理(CRT)(壮哉我大中华╰(*°▽°*)╯)
niiick的中国剩余定理 && 扩展中国剩余定理
并不是很难的东西,原理跟多元线性同余方程类似,扩展中国剩余定理更是同一个东西,只不过推导的过程有些不一样。这里便不重复。
BiorhythmsPOJ - 1006
P3868 [TJOI2009]猜数字
P4777 【模板】扩展中国剩余定理(EXCRT)
Strange Way to Express Integers POJ - 2891
#include<cstdio> typedef long long ll; ll bb[5],cc[5]={23,28,33},cp; ll exgcd(ll a,ll b, ll &x,ll &y){ if(!b){ x=1; y=0; return a; } ll g=exgcd(b,a%b,y,x); y-=a/b*x; return g; } ll inv(ll a,ll c){ ll g,x,y; g=exgcd(a,c,x,y); return g==1 ? (x%c+c)%c : -1; } ll crt(int n){ ll ans=0,temp; cp=1; for(int i=0;i<n;i++) cp*=cc[i]; for(int i=0;i<n;i++){ temp=cp/cc[i]; ans+=bb[i]*temp*inv(temp,cc[i]); if(ans>=cp) ans%=cp; } return (ans+cp)%cp; } int main(){ int t=1; ll d,ans; while(~scanf("%lld%lld%lld%lld",&bb[0],&bb[1],&bb[2],&d)){ if(bb[0]==-1&&bb[1]==-1&&bb[2]==-1) break; ans=crt(3)-d; if(ans<=0) ans+=cp; printf("Case %d: the next triple peak occurs in %lld days. ",t++,ans); } return 0; }
1 #include<cstdio> 2 typedef long long ll; 3 const int N=1e5+11; 4 ll bb[N],cc[N]; 5 ll exgcd(ll a,ll b, ll &x,ll &y){ 6 if(!b){ 7 x=1; 8 y=0; 9 return a; 10 } 11 ll g=exgcd(b,a%b,y,x); 12 y-=a/b*x; 13 return g; 14 } 15 ll excrt(int n){ 16 ll ans=0,cp=1,a,b,c,g,x,y; 17 for(int i=0;i<n;i++){ 18 a=cp,c=cc[i],b=bb[i]-ans; 19 g=exgcd(a,c,x,y); 20 if(b%g) return -1; 21 b/=g;c/=g; 22 x=(x*b%c+c)%c; 23 ans+=x*cp; 24 cp*=c; 25 ans=(ans%cp+cp)%cp; 26 } 27 return ans; 28 } 29 int main(){ 30 int n; 31 while(~scanf("%d",&n)){ 32 for(int i=0;i<n;i++) scanf("%lld%lld",&cc[i],&bb[i]); 33 printf("%lld ",excrt(n)); 34 } 35 return 0; 36 }