转载--数论总结

一.素数打表：

1.平常打表：

[cpp]  view plain copy

1. for(i=2;i<=n;i++)      
2.         if(!s[i])      
3.         {      
4.             for(j=2*i;j<=n;j+=i)      
5.                 s[j]=1;      
6.         }   

2.线性打表：

[cpp]  view plain copy

1. void get_prime()    
2. {    
3.     int cnt = 0;    
4.     for (int i = 2; i < N; i++)    
5.     {    
6.         if (!tag[i])    
7.             p[cnt++] = i;    
8.         for (int j = 0; j < cnt && p[j] * i < N; j++)    
9.         {    
10.             tag[i*p[j]] = 1;    
11.             if (i % p[j] == 0)    
12.                 break;    
13.         }    
14.     }    
15. }    

二.欧几里得与扩展欧几里得算法：

1.欧几里得：

[cpp]  view plain copy

1. int gcd(int a,int b)  
2. {  
3.     return b==0?a:gcd(b,a%b);  
4. }  

2.扩展欧几里得：

[cpp]  view plain copy

1. //求整数x和y，使得ax+by=d，且|x|+|y|最小。其中d=gcd(a,b)  
2. void extend_gcd(int a,int b,int &d,ll &x,ll &y)  
3. {  
4.     if(!b) d=a,x=1,y=0;  
5.     else {extend_gcd(b,a%b,d,y,x);y-=x*(a/b);}  
6. }  

三.快速幂取模：

[cpp]  view plain copy

1. //求a^b%n  
2. ll mod(ll a,ll b,int n)  
3. {  
4.     if(b==0) return 1;  
5.     ll ans=mod(a,b/2,n);  
6.     ans=ans*ans%n;  
7.     if(b&1) ans=ans*a%n;  
8.     return ans;  
9. }  

四.欧拉phi函数：

欧拉函数phi(n)等于不超过n且和n互素的整数个数。 互素：两个数的最大公约数为1的数称为互素数。

1.单个欧拉函数求法代码：

[cpp]  view plain copy

1. int euler_phi(int n)  
2. {  
3.     int m=sqrt(n+0.5),ans=n;  
4.     for(int i=2;i<=m;i++)  
5.         if(n%i==0)  
6.         {  
7.             ans=ans/i*(i-1);  
8.             while(n%i==0) n/=i;  
9.         }  
10.     if(n>1) ans=ans/n*(n-1);  
11.     return ans;  
12. }  

2.函数表：多个欧拉函数一起求，类似筛选法计算phi(1),phi(2),……phi(n).

[cpp]  view plain copy

1. int phi[MM]  
2. void phi_table(int n)  
3. {  
4.     mem(phi,0); phi[1]=1;  
5.     for(int i=2;i<=n;i++)  
6.         if(!phi[i])  
7.             for(int j=i;j<=n;j+=i)  
8.             {  
9.                 if(!phi[j]) phi[j]=j;  
10.                 phi[j]=phi[j]/i*(i-1);  
11.             }  
12. }  

五.乘法逆：

在n的同余类中的两个数a和b满足ab=1,比如在15的同余类中7*13=1，在这种情况下，说a和b互为乘法的逆。即7的逆的13，13的逆为7；如果知道a就可以求b，知道b也就可以求a。代码如下：

1.用扩展欧几里得求：

[cpp]  view plain copy

1. //计算模n下a的逆。如果不存在逆，返回-1  
2. int inv(int a,int n)  
3. {  
4.     int d,x,y;  
5.     extend_gcd(a,n,d,x,y); //扩展欧几里得函数  
6.     return d==1?(x+n)%n:-1;  
7. }  

2.利用欧拉定理求：

[cpp]  view plain copy

1. //a的逆就是mod(a,n-2,n)  
2. ll mod(ll a,ll n-2,int n)  
3. {  
4.     if(b==0) return 1;  
5.     ll ans=mod(a,b/2,n);  
6.     ans=ans*ans%n;  
7.     if(b&1) ans=ans*a%n;  
8.     return ans;  
9. }  

六.模方程：

1.线性模方程：axºb(mod n)

给a线性模方程：ax ºb(mod n)。求出x,因为x有多个解，所以代码也有多个解。把其化成ax-ny=b，d=gcd(a,n)不是b的约数时无解，否则两边同时除以d……即利用扩展欧几里得求解。

[cpp]  view plain copy

1. //返回0时不存在解，为1有解  
2. int mod_gcd(int a,int b,int n)  
3. {  
4.     int x,y,d,x0,i;  
5.     extend_gcd(a,n,d,x,y);  
6.     if(b%d) return 0;  
7.     x0=x*(b/d)%n;  
8.     for(i=1;i<d;i++)  
9.         printf("%d\n",(x0+i*(n/d))%n);  
10.     return 1;  
11. }  

2.费马小定理：a^(p-1) ≡1（mod p）

费马小定理是数论中的一个重要定理，其内容为： 假如p是质数，且(a,p)=1，那么 a^(p-1) ≡1（mod p）。

即：假如p是质数，且a,p互质，那么a的(p-1)次方除以p的余数恒等于1。

3.中国剩余定理：xºai(mod mi)

如果线性模方程有多个，xºai(mod mi)该怎么做呢？这里就可以用中国剩余定理了。令M为所有mi的乘积，则在模M的剩余系下原方程组有唯一解。

但是用中国剩余定理的条件是mi之间是两两互质的数：

[cpp]  view plain copy

1. //n个方程：x=a[i](mod m[i]) (0<=i<n)  
2. ll china(int n,int *a,int *m)  
3. {  
4.     ll M=1,d,y,x=0;  
5.     for(int i=0;i<n;i++)  
6.         M*=m[i];  
7.     for(int i=0;i<n;i++)  
8.     {  
9.         ll w=M/m[i];  
10.         gcd(m[i],w,d,d,y);  
11.         x=(x+y*w*a[i])%M;  
12.     }  
13.     return (x+M)%M;  
14. }  

如果mi不是两两互质的数，那该怎么办呢？那就得要扩展欧几里得算法把两个等式合并成一个了，然后继续下去求了……（POJ 2891)

[cpp]  view plain copy

1. ll gcd(ll a,ll b)    
2. {    
3.     ll t=(!b?a:gcd(b,a%b));   //记忆递归省点时间    
4.     return t;    
5. }    
6. void exgcd(ll a,ll b,ll &x,ll &y)    
7. {    
8.     if(!b) {x=1;y=0;return;}    
9.     exgcd(b,a%b,y,x);  //这里用记忆化递归时间还是一样的，所以不用    
10.     y-=a/b*x;    
11. }    
12. int main()    
13. {    
14.     ll i,k,a1,r1,a2,r2,c,x,y,l,t;    
15.     while(cin>>k)    
16.     {    
17.         int flag=0;    
18.         cin>>a1>>r1;    
19.         for(i=1;i<k;i++)    
20.         {    
21.             cin>>a2>>r2;    
22.             if(flag) continue;    
23.             c=r2-r1;    
24.             l=gcd(a1,a2);    
25.             if(c%l!=0) {flag=1;continue;}    
26.             exgcd(a1,a2,x,y);    
27.             t=a2/l;    
28.             x=((c/l*x)%t+t)%t;    
29.             r1+=a1*x;    
30.             a1*=t;    
31.         }    
32.         if(flag) cout<<-1<<endl;    
33.         else cout<<r1<<endl;    
34.     }    
35.     return 0;    
36. }    

4.高次同余模方程求解：a^x º b(mod n)

利用欧拉定理求解模方程，当n为素数时，解模方程a^xºb(mod n)。根据欧拉定理，只需检查x=0,1,2,……,n-1是不是解即可因为a^(n-1)º1(mod n),当x超过n-1时a^x就开始循环了。

[cpp]  view plain copy

1. //求解模方程a^xb(mod n)。n为素数。无解返回0  
2. int log_mod(int a,int b,int n)  
3. {  
4.     int m,v,e=1,i;  
5.     m=sqrt(n+0.5);  
6.     v=inv(mod(a,m,n),n);  
7.     map<int,int>x;  
8.     x[1]=0;  
9.     for(i=1;i<m;i++)  
10.     {  
11.         e=e*a%n;  
12.         if(!x.count(e)) x[e]=i;  
13.     }  
14.     for(i=0;i<m;i++)  
15.     {  
16.         if(x.count(b)) return i*m+x[b];  
17.         b=b*v%n;  
18.     }  
19.     return 0;  
20. }  

哈希法的高次同余方程模板：优点快速，代码少。

[cpp]  view plain copy

1. #define MOD 76543  
2. int hs[MOD],head[MOD],next[MOD],id[MOD],top;  
3. void insert(int x,int y)  
4. {  
5.     int k = x%MOD;  
6.     hs[top] = x, id[top] = y, next[top] = head[k], head[k] = top++;  
7. }  
8. int find(int x)  
9. {  
10.     int k = x%MOD;  
11.     for(int i = head[k]; i != -1; i = next[i])  
12.         if(hs[i] == x)  
13.             return id[i];  
14.     return -1;  
15. }  
16. int BSGS(int a,int b,int n) //a^L=b(mod n) 求L  
17. {  
18.     memset(head,-1,sizeof(head));  
19.     top = 1;  
20.     if(b == 1)return 0;  
21.     int m = sqrt(n*1.0), j;  
22.     long long x = 1, p = 1;  
23.     for(int i = 0; i < m; ++i, p = p*a%n)insert(p*b%n,i);  
24.     for(long long i = m; ; i += m)  
25.     {  
26.         if( (j = find(x = x*p%n)) != -1 )return i-j;  
27.         if(i > n)break;  
28.     }  
29.     return -1;  //L不存在的情况  
30. }  

还有就是算法经典入门里的模板，用的是map<int,int>x，缺点：太慢，而且代码又多！在POJ 2417上TLE了。

[cpp]  view plain copy

1. void extend_gcd(ll a,ll b,ll &d,ll &x,ll &y)  
2. {  
3.     if(!b) d=a,x=1,y=0;  
4.     else {extend_gcd(b,a%b,d,y,x);y-=x*(a/b);}  
5. }  
6. ll mod(ll a,ll b,ll n)  
7. {  
8.     if(b==0) return 1;  
9.     ll ans=mod(a,b/2,n);  
10.     ans=ans*ans%n;  
11.     if(b&1) ans=ans*a%n;  
12.     return ans;  
13. }  
14. ll inv(ll a,ll n)  
15. {  
16.     ll d,x,y;  
17.     extend_gcd(a,n,d,x,y); //扩展欧几里得函数  
18.     return d==1?(x+n)%n:-1;  
19. }  
20. ll log_mod(ll a,ll b,ll n)  
21. {  
22.     ll m,v,e=1,i;  
23.     m=sqrt(n+0.5);  
24.     v=inv(mod(a,m,n),n);  
25.     map<int,int>x;  
26.     x[1]=0;  
27.     for(i=1;i<m;i++)  
28.     {  
29.         e=e*a%n;  
30.         if(!x.count(e)) x[e]=i;  
31.     }  
32.     for(i=0;i<m;i++)  
33.     {  
34.         if(x.count(b)) return i*m+x[b];  
35.         b=b*v%n;  
36.     }  
37.     return -1;  
38. }  

七.Lucas定理：

A、B是非负整数，p是质数。A B写成p进制：A=a[n]a[n-1]...a[0]，B=b[n]b[n-1]...b[0]。
则组合数C(A,B)与C(a[n],b[n])*C(a[n-1],b[n-1])*...*C(a[0],b[0])  mod p同余

即：Lucas(n,m,p)=C(n%p,m%p)*Lucas(n/p,m/p,p)

求大组合数取模：C(n,m)%p，大组合数取模。

[cpp]  view plain copy

1. ll fac[100003];  
2. void init(ll p)  
3. {  
4.     fac[0] = 1;  
5.     for (int i=1; i<=p; i++) fac[i] = fac[i-1]*i%p;  
6. }  
7. ll PowerMod(ll a, ll b, ll k)  
8. {  
9.     ll tmp = a, ret = 1;  
10.     while (b)  
11.     {  
12.         if (b & 1) ret = ret * tmp % k;  
13.         tmp = tmp * tmp % k;  
14.         b >>= 1;  
15.     }  
16.     return ret;  
17. }  
18. ll Lucas(ll n, ll m, ll p)  
19. {  
20.     ll ret = 1;  
21.     while (n && m)  
22.     {  
23.         ll nn = n%p, mm = m%p;  
24.         if (nn < mm) return 0;  
25.         ret = ret*fac[nn]*PowerMod(fac[mm]*fac[nn-mm]%p, p-2, p)%p;  
26.         n /= p;  
27.         m /= p;  
28.     }  
29.     return ret;  
30. }  

八.指数循环节：

即求：a[0]^a[1]^……^a[n]的值，比如：2^4^3^5=2。

[cpp]  view plain copy

1. //求a[0]^a[1]^……^a[n]%m的值（指数循环节）  
2. LL a[20],n,m;  
3. LL euler(LL n)  
4. {  
5.     LL m=(LL)sqrt(n+0.5);  
6.     LL ans=n;  
7.     for(LL i=2; i<=m; i++)  
8.         if(n%i==0)  
9.         {  
10.             ans=ans/i*(i-1);  
11.             while(n%i==0)n/=i;  
12.         }  
13.     if(n>1)ans=ans/n*(n-1);  
14.     return ans;  
15. }  
16. LL pmod(LL a,LL n,LL m)  
17. {  
18.     LL now = 1;  
19.     for(LL i=0; i<n; i++)  
20.     {  
21.         now *= a;  
22.         if(now>=m)  break;  
23.     }  
24.     if(now >= m)   now = m;  
25.     else now=0;  
26.   
27.     if(n==0)return 1;  
28.     LL x=pmod(a,n/2,m);  
29.     LL ans=(LL)x*x%m;  
30.     if(n%2==1) ans=ans*a%m;  
31.     return ans+now;  
32. }  
33. LL solve(LL cur,LL n,LL m)  
34. {  
35.     if(cur==n-1)  
36.     {  
37.         if(a[cur]>=m)  
38.             return a[cur]%m+m;  
39.         else  
40.             return a[cur];  
41.     }  
42.     LL M=euler(m);  
43.     LL p=solve(cur+1,n,M);  
44.     LL ans=pmod(a[cur],p,m);  
45.     return ans;  
46. }  
47. //输出的是solve(0,n,m)%m  

九.高斯消元法：

[cpp]  view plain copy

1. //要求系数矩阵可逆  
2. //这里的A是增广矩阵，即A[i][n]是第i个方程右边的常数bi。  
3. //运行结束后A[i][n]是第i个未知数的值  
4. typedef double Matrix[MM][MM];  
5. void gauss(Matrix A,int n)  
6. {  
7.     int i,j,k,r;  
8.     //消元过程  
9.     for(i=0;i<n;i++)  
10.     {  
11.         //选一行r并与第i交换  
12.         r=i;  
13.         for(j=i+1;j<n;j++)  
14.             if(fabs(A[i][j])>fabs(A[r][i])) r=j;  
15.         if(r!=i)  
16.         {  
17.             for(j=0;j<=n;j++)  
18.                 swap(A[r][j],A[i][j]);  
19.         }  
20.         //与第i+1~n行进行消元  
21.         for(j=n;j>=i;j--)  
22.             for(k=i+1;k<n;k++)  
23.                 A[k][j]-=A[k][i]/A[i][i]*A[i][j];  
24.     }  
25.     for(i=n-1;i>=0;i--)  
26.     {  
27.         for(j=i+1;j<n;j++)  
28.             A[i][n]-=A[j][n]*A[i][j];  
29.         A[i][n]/=A[i][i];  
30.     }  
31. }  

十.大素数判断Miller_Rabin 算法：

判断的句子是：if(Is_Prime(n)) cout<<"YES"<<endl;

[cpp]  view plain copy

1. const int TIMES=200;  //判断次数可以减少，但是越多，判断数组的准确率越高  
2. long long Gcd(long long a,long long b)  
3. {  
4.     if(b==0)  
5.     {  
6.         return a;  
7.     }  
8.     else  
9.     {  
10.         return Gcd(b,a%b);  
11.     }  
12. }  
13. long long Cheng_Mod(long long a,long long x,long long p)  
14. {  
15.     long long Ret=0,s=a;  
16.     while(x>0)  
17.     {  
18.         if(x%2==1)  
19.         {  
20.             Ret=(Ret+s)%p;  
21.         }  
22.         s=(s+s)%p;  
23.         x/=2;  
24.     }  
25.     return Ret;  
26. }  
27. long long Mult_Mod(long long a,long long x,long long p)  
28. {  
29.     long long Ret=1,s=a;  
30.     while(x>0)  
31.     {  
32.         if(x%2==1)  
33.         {  
34.             Ret=Cheng_Mod(Ret,s,p);  
35.         }  
36.         s=Cheng_Mod(s,s,p);  
37.         x/=2;  
38.     }  
39.     return Ret;  
40. }  
41. bool Is_Prime(long long p)  
42. {  
43.     if(p==1)  
44.     {  
45.         return 0;  
46.     }  
47.     if(p<=3)  
48.     {  
49.         return 1;  
50.     }  
51.     int i,j;  
52.     long long b=0,m=p-1,d=p-1,x;  
53.     while(m%2==0)  
54.     {  
55.         b++;  
56.         m/=2;  
57.     }  
58.     for(i=1; i<=TIMES; i++)  
59.     {  
60.         x=rand()%(p-3)+2;  
61.         d=Mult_Mod(x,m,p);  
62.         if(d==1||d==p-1)  
63.         {  
64.             continue;  
65.         }  
66.         for(j=1; j<=b; j++)  
67.         {  
68.             d=Cheng_Mod(d,d,p);  
69.             if(d==1)  
70.             {  
71.                 return 0;  
72.             }  
73.             if(d==p-1)  
74.             {  
75.                 break;  
76.             }  
77.         }  
78.         if(d!=p-1)  
79.         {  
80.             return 0;  
81.         }  
82.     }  
83.     return 1;  
84. }