数论知识总结——史诗大作(这是一个flag)

1、快速幂

计算a^b的快速算法,例如,3^5,我们把5写成二进制101,3^5=3^1*1+3^2*2+3^4*1

1 ll fast(ll a,ll b){ll ans=1;for(;b;b>>=1,a=mul(a,a))if(b&1)ans=mul(ans,a);return ans;}//一行快速幂
View Code

 

 

2、快速乘

当模数较大时,直接乘会爆掉long long,需要快速乘法。

即用浮点计算倍数,做差相当于计算余数模2^63的结果,然后再模一下就好了(因为余数不超过long long)

1 typedef long long ll;
2 ll mul(ll x,ll y){return ((x*y-(ll)(((long double)x*y+0.5)/mod)*mod)%mod+mod)%mod;}//一行快速乘
View Code

练习题目:【HDU5187】zhx's contest

 

 

3、同余原理(gcd)

定理:gcd(a,b)=gcd(b,a mod b)

证明:a可以表示成a=kb+r,则r=a mod b 

     假设d是a,b的一个公约数,则有 d|a, d|b,而r = a - kb,因此d|r ,因此d是(b,a mod b)的公约数

     假设d 是(b,a mod b)的公约数,则 d|b , d|r ,而a=kb+r,因此d也是(a,b)的公约数

     因此(a,b)和(b,a mod b)的公约数是一样的,其最大公约数也必然相等,得证。

1 int gcd(a,b){return !b?a:gcd(b,a%b);}//一行gcd
View Code

 

 

4、丢番图方程

裴蜀定理:丢番图方程(二元一次,下同)ax+by=m有解当且仅当m|(a,b)

扩展gcd:求丢番图方程ax+by=gcd(a,b)的整数解

证明:a*x1+b*y1=gcd(a,b)

   b*x2+(a mod b)*y2=gcd(b,a mod b)

     因为gcd(a,b)=gcd(b,a mod b)

   得a*x1+b*y1= b*x2+(a mod b)*y2 = b*x2+(a-a/b*b)*y2 = a*y2+b*(x2-a/b*y2)

   所以x1=y2,y1=x2-a/b*y2

末状态:b=0,a=gcd(a,b)时,gcd(a,b)*x=gcd(a,b),得x=1

扩展欧几里得的过程可以理解为从末状态向上不断回溯的过程,直到得到原方程的一组解。

1 void exgcd(int a,int b,int &x,int &y)
2 {
3     if(b==0)  {x=1; y=0; return;}
4     exgcd(b,a%b,x,y);
5     int t=x;x=y;y=t-a/b*y;
6 }
View Code

丢番图方程的通解:若(a,b)=d, ax+by=m有一组解(x0,y0)通解为:x=x0+(b/d)k, y=y0-(a/d)k就是设法让正负相互对消,(真·小学奥数内容)

解惑:多数初学者可能会有这样的疑问(反正我刚学的时候有),要解的方程是ax+by=m,而以上算法解的是ax+by=gcd(a,b)

其实ax+by=gcd(a,b)可以变为ax*k+by*k=gcd(a,b)*k,令gcd(a,b)*k=m,求出k,然后x*k就是ax+by=m的解。

练习题目:【NOIP2012】同余方程

【POJ1061】青蛙的约会

【NOI2002】荒岛野人

【CF#724C】Ray Tracing

 

 

5、乘法逆元

定义:对于正整数a和m,如果a*x mod m=1, 则x的最小正整数解称为a模m的逆元,表示为a^(-1) (mod m)

求解方法:(1)根据费马小定理:有a^(m-1) mod m=1,  即a*a^(m-2) mod m=1, 所以a^(m-2)就是a模m的逆元。 (适用条件:m为素数)

     (2)扩展欧几里得求解:a*x mod m=1,即a*x-m*y=1,解这个丢番图方程即可。  (适用条件:gcd(a,m)=1)

     (3)欧拉函数法:留个坑。

     (4)线性递推:规定m为质数,且1^(-1)≡1(mod m)    设m=k∗a+b  (b<a,1<a<m),  即 k∗a+b≡0(mod m) 

              两边同时乘以a^(−1)∗b^(−1),得到 k∗b^(−1)+a^(−1)≡0   (mod m) 

              a^(-1)≡−k∗b^(−1)     (mod m)      

              a^(−1)≡−m/a∗(m mod a)^(-1)      (mod m) 

              从头开始扫一遍即可,时间复杂度O(n)         (适用条件:m为素数)

逆元的一个经典应用:我们知道除法不符合模性质,那么在解决(a/b)mod m 的问题时该怎么办呢?

     推导一下:设a/b=k*m+x,那么x就是所求答案。

          a=k*m*b+x*b

          a mod (m*b) = x*b

          x=(a mod (m*b))/b

但是如果b很大,则会出现爆精度问题,所以我们避免使用除法直接计算。 

可以使用逆元将除法转换为乘法: 

假设b存在乘法逆元,即与m互质(充要条件)。设c是b的逆元,即b∗c≡1(modm),那么有a/b=(a/b)∗1=((a/b)∗b∗c)mod m=(a∗c)mod m

即,除以一个数取模等于乘以这个数的逆元取模。

 

 

6、中国剩余定理

<1>互质版

定理内容:(我不会证明)

给定n组同余方程(m1,m2,m3……mn两两互质):

x≡a1(mod m1)

x≡a2(mod m2)

x≡a3(mod m3)

……

x≡an(mod mn)

令M=m1∗m2∗...∗mn , 则x=(a1*M1*M^(-1)+a2*M2*M^(−1)+...+an*Mn*M^(−1)) mod M (其中Mi=M/mi , M^(-1)为Mi mod mi的乘法逆元。)

1 void China()
2 {
3     for(ll i=1;i<=n;i++)
4     {
5         ll x,y,Mi=M/m[i];
6         exgcd(Mi,m[i],x,y);
7         ans=(ans+Mi*x*a[i])%M;
8     }
9 }
View Code

<2>非互质版

设x%a1=m1, x%a2=m2

则x=k1*m1+a1=k2*m2+a2

设(m1,m2)=d,则k1*(m1/d) = k2*(m2/d)+(a2-a1)/d   (注意必须满足d|(a2-a1))

注意到m1/d与m2/d互质,解方程:k1*(m1/d)=(a2-a1)/d (mod m2/d)得到通解k1 = k+r*(m2/d)

代回得x = (k+r*m2/d)*m1+a1 = r*(m1*m2/d)+k*m1+a1

即x=k*m1+a1 (mod lcm[m1,m2])

这样就实现了2个方程的合并,将所有方程两两合并,就能求出最终的x

 1 int China()
 2 {
 3     int A=a[1],k,y;  M=m[1];
 4     for(int i=2;i<=n;i++)
 5     {
 6         int da=a[i]-A,g;
 7         exgcd(M,m[i],g,k,y);
 8         if(da%g)  return -1; 
 9         k*=da/g;
10         A+=k*M;
11         M=M*m[i]/g;
12         A=(A+M)%M;
13     }
14     return A;
15 }
View Code

练习题目:【Vijos1164】曹冲养猪

【codevs3990】中国余数定理2

【HDU1573】X问题

【POJ2891】Strange Way to Express Integers

【CF#338D】GCD Table

 

 

7、欧拉函数

定义:对于一个正整数n,小于n且和n互质的正整数(包括1)的个数叫做欧拉函数,记作φ(n) 。

通式:φ(x)=x*(1-1/p1)*(1-1/p2)*(1-1/p3)*(1-1/p4)…..(1-1/pn),其中p1, p2……pn为x的所有质因数,x是不为0的整数。φ(1)=1(唯一和1互质的数就是1本身)。

欧拉函数是积性函数——若m,n互质,φ(mn)=φ(m)φ(n)。

特殊性质:1、当n为奇数时,φ(2n)=φ(n)

     2、M=d|Mϕ(d) (d|M)

     3、对于给定的一个素数p,φ(p)=p -1。则对于正整数 n = p^k ,φ(n) = p^k - p^(k -1)

递推性质:若(N % a == 0 && (N / a) % a == 0) 则有:E(N)=E(N / a) * a;

       若(N % a == 0 && (N / a) % a != 0)  则有:E(N)=E(N / a) * (a - 1)。

由递推性质可得求欧拉函数的代码,时间复杂度:O(n)

 1 void get()
 2 {
 3     memset(check,0,sizeof(check));
 4     int cnt=0;
 5     for(int i=2;i<=N;i++)
 6     {
 7         if(!check[i])  {prime[++cnt]=i;  phi[i]=i-1;}
 8         for(int j=1;j<=cnt&&prime[j]*i<=N;j++)
 9         {
10             check[i*prime[j]]=1;
11             if(i%prime[j])  phi[i*prime[j]]=phi[i]*(prime[j]-1);
12             else  {phi[i*prime[j]]=phi[i]*prime[j]; break;}
13         }
14     }
15 }
View Code

练习题目:【bzoj2190】[SDOI2008]仪仗队          

【POJ2478】Farey Sequence                

【bzoj2818】Gcd               

【bzoj2705】[SDOI2012]Longge的问题     

【bzoj2186】[Sdoi2008]沙拉公主的困惑    

【bzoj2749】[HAOI2012]外星人

 

 

8、数论四大定理

<1>欧拉定理

内容:对于互质的正整数a和n,有a*φ(n) ≡ 1 (mod n)。

用欧拉定理解决幂的模运算:a^b mod m = a^(b mod phi[m]) mod m

证明:a^b mod m=a^((phi[m])*k)*a^(b mod phi[m]) mod m=a^(b mod phi[m]) mod m

练习题目:【bzoj3884】上帝与集合的正确用法

<2> 费马小定理

内容:a^(p-1)%p=1  (p是素数)

证明:根据欧拉定理有a^(phi[p])%p=1, 因为p是质数,所以phi[p]=p-1,所以a^(p-1)%p=1 

<3>威尔逊定理

内容: (p-1)! mod p=-1 (p是素数)

证明:考虑一个长度为p的环,将环上p个点顺次连成一条折线

     共有(p-1)!种,其中p-1种全同的具有循环不变性,其他每个都属于大小为p的同余类

     故p|((p-1)!-(p-1))

     即(p-1)!=-1(mod p)

练习题目:【HUD5391】Zball in Tina Town

<4>中国剩余定理

上面讲过了,不再赘述。

 

 

9、素数相关

<1>线性筛素数

可以在O(n)的时间内求出1~n的所有素数。

 1 int cnt,prime[MAXN],check[MAXN];
 2 void get()
 3 {
 4     for(int i=2;i<=n;i++)
 5     {
 6         if(!check[i])  prime[++cnt]=i;
 7         for(int j=1;j<=cnt&&prime[j]*i<=n;j++)
 8         {
 9             check[prime[j]*i]=1;
10             if(i%prime[j]==0)  break;//保证线筛的时间复杂度
11         }
12     }
13 }
View Code

<2>Rabin-Miller测试

根据费马小定理,若p为素数,则必有a^(p-1) mod p=1 对和p互质的a成立。

二次探测定理:如果p是素数,且0<x<p,则方程x^2 mod p=1的解为1或p-1。

    证明:因为x^2 mod p=1,所以p|(x-1)*(x+1),故(1%p)的平方根为1或-1

所以若p为素数,则必有a^(p-1) mod p 的平方根为1或-1

分解p-1为d*2^s,其中d为奇数

从i=0逐次计算a^(d*2^(s-i)),相当于“开方”,若得到-1或追查到a^d=1 (mod p),则p通过测试,否则不通过

时间复杂度O(k*(logn)^3),k是选的a的个数

 1 int fast(int a,int b,int mod) {int sum=1;for(;b;b>>=1,a=a*a%mod)if(b&1)sum=sum*a%mod;return sum;}//一行快速幂
 2 bool Rabin_Miller(int p,int a)
 3 {
 4     if(p==2)  return 1;
 5     if(p&1==0||p==1)  return 0;
 6     int d=p-1;
 7     while(!(d&1))  d>>=1;
 8     int m=fast(a,d,p);
 9     if(m==1)  return 1;
10     for(;d<p;d<<=1,m=m*m%p) if(m==p-1) return 1;
11     return 0;
12 }
13 bool isprime(int x)
14 {
15     static int prime[9]={2,3,5,7,11,13,17,19,23};//相当于全局变量
16     for(int i=0;i<9;i++)
17     {
18         if(x==prime[i])  return 1;
19         if(!Rabin_Miller(x,prime[i]))  return 0;
20     }
21     return 1;
22 }
View Code

练习题目:【HUD2138】How many prime numbers

<3>Pollard_Rho因数分解 

基础:一个能生成随机序列的函数f,一般取f(x)=x^2+c  (c是一个随机的数)

递归处理:

若n是素数(用Rabin-Miller判断),则已完成分解

否则,从某个x1开始,不断计算x1=f(x1)%n

在模n的任一因数意义下,它们将ρ型循环

有一堆ρ型,故名Pollard-Rho

若出现循环,即1<gcd(abs(x1-x2),n)<n,则找到了因数gcd(abs(x1-x2),n)

递归分解gcd(abs(x1-x2),n)和n/gcd(abs(x1-x2),n)

若出现重复,则更改函数(更改c的值)

判重复:判断x1和x[2..2]是否检出,x2和x[3..4]是否检出,x4和x[5..8]是否检出……

时间复杂度:找出因子p需要O(sqrt(p))

 1 int cnt,fac[1000100],num[1000100];
 2 int gcd(int a,int b) {return !b?a:gcd(b,a%b);}
 3 int fast(int a,int b,int mod) {int sum=1;for(;b;b>>=1,a=a*a%mod)if(b&1)sum=sum*a%mod;return sum;}//一行快速幂
 4 bool Rabin_Miller(int p,int a)
 5 {
 6     if(p==2)  return 1;
 7     if(p&1==0||p==1)  return 0;
 8     int d=p-1;
 9     while(!(d&1))  d>>=1;
10     int m=fast(a,d,p);
11     if(m==1)  return 1;
12     for(;d<p;d<<=1,m=m*m%p) if(m==p-1) return 1;
13     return 0;
14 }
15 bool isprime(int x)
16 {
17     static int prime[9]={2,3,5,7,11,13,17,19,23};//相当于全局变量
18     for(int i=0;i<9;i++)
19     {
20         if(x==prime[i])  return 1;
21         if(!Rabin_Miller(x,prime[i]))  return 0;
22     }
23     return 1;
24 }
25 void Pollord_Rho(int x)
26 {
27     if(isprime(x))  {fac[++cnt]=x;  return;}
28     int c=3;
29     while(1)
30     {
31         int x1(1),x2(1),k(2),i(1);
32         while(1)
33         {
34             x1=((x1*x1)%x+c)%x;
35             int d=gcd(abs(x1-x2),x);
36             if(d>1&&d<x)
37             {
38                 Pollord_Rho(d);
39                 Pollord_Rho(x/d);
40                 return;
41             }
42             if(x1==x2)  break;
43             if(++i==k)  k<<=1,x2=x1;
44         }
45         c++;
46     }
47 }
48 void solve(int x)
49 {
50     Pollord_Rho(x);
51     sort(fac+1,fac+cnt+1);
52     int len=0;
53     for(int i=1;i<=cnt;i++)
54     {
55         if(fac[i]==fac[i-1])  num[len]++;
56         else num[++len]=1,fac[len]=fac[i];
57     }
58     for(int i=1;i<len;i++)  printf("%d^%d*",fac[i],num[i]);
59     printf("%d^%d\n",fac[len],num[len]);
60 }
View Code

练习题目:【bzoj3667】Rabin-Miller算法 

【POJ1811】Prime Test

 

 

10、离散对数

离散对数是解决同余方程 A^x = B (mod C)的最小整数解

<1>Shank大步小步法(又叫Baby Step Giant Step算法

由欧拉定理,只需要检查1<=x<phi(M),检查到M也行  

的值存到一个Hash表里面

的值一一枚举出来,每枚举一个就在Hash表里面寻找是否有一个值满足

,如果有则找到答案,否则继续

最终答案就是的值对应的原来A的幂

那么如果C不是素数呢?

当C为合数时,我们就需要把Baby Step Giant Step扩展一下。在普通Baby Step Giant Step中,由于是素数,那么,所以一定有唯一解的。那么,当为合数时,我们可以这样处理:

对于方程,我们拿出若干个A出来与C来消去公共因子,使得为止,那么此时我们就可以直接通过扩展欧几里得来计算结果了。

——以上内容转自ACdreamers离散对数(Baby Step Giant Step)

下面给出bzoj2480的代码作为模板(良心的博主加上了注释):

 1 /*************
 2   bzoj 2480
 3   by chty
 4   2016.11.8
 5 *************/
 6 #include<iostream>
 7 #include<cstdio>
 8 #include<cstring>
 9 #include<cstdlib>
10 #include<ctime>
11 #include<cmath>
12 #include<algorithm>
13 using namespace std;
14 #define mod 99991
15 typedef long long ll;
16 struct node{ll v,num,f;}hash[mod+10];
17 ll A,B,C;
18 inline ll read()
19 {
20     ll x=0,f=1;  char ch=getchar();
21     while(!isdigit(ch))  {if(ch=='-')  f=-1;  ch=getchar();}
22     while(isdigit(ch))  {x=x*10+ch-'0';  ch=getchar();}
23     return x*f;
24 }
25 ll gcd(ll a,ll b) {return !b?a:gcd(b,a%b);}
26 void insert(ll v,ll x)  //哈希表插入操作
27 {
28     ll t=v%mod;
29     while(hash[t].f&&hash[t].v!=v) {t++; if(t>mod) t-=mod;}
30     if(!hash[t].f) {hash[t].f=1;hash[t].num=x;hash[t].v=v;}
31 }
32 ll find(ll v)  //哈希表查询操作
33 {
34     ll t=v%mod;
35     while(hash[t].f&&hash[t].v!=v) {t++; if(t>mod) t-=mod;}
36     if(!hash[t].f)  return -1;
37     else return hash[t].num;
38 }
39 void exgcd(ll a,ll b,ll &x,ll &y)  //扩展欧几里得算法求丢番图方程的一组解
40 {
41     if(!b)  {x=1; y=0; return;}
42     exgcd(b,a%b,x,y);
43     ll t=x;x=y;y=t-a/b*y;
44 }
45 ll Shank()
46 {
47     ll ret=1;
48     for(ll i=0;i<=50;i++) {if(ret==B)return i;ret=ret*A%C;}  
49     //判断在0~50内是否有解,这就是Baby_Step,不加这一步会wa掉,因为下面做除法时,假设的是A足够用
50     ll temp,ans(1),cnt(0);
51     while((temp=gcd(A,C))!=1)//当A与C不互质时,拿出若干个A出来与C来消去公共因子,使得gcd(A,C)=1为止
52     {
53         if(B%temp)  return -1;//根据裴蜀定理,如果gcd(A,C)|B不成立,则无解
54         C/=temp;  B/=temp;
55         ans=ans*(A/temp)%C;  
56         cnt++;  //记录已经乘了cnt个A,以及这些A与C消掉公因子后的结果ans
57     }
58     ll m=(ll)ceil(sqrt(C*1.0)),t(1);  //保证时间复杂度为O(n^0.5)
59     for(ll i=0;i<m;i++) {insert(t,i);t=t*A%C;}  //把A^0~A^m存到哈希表中
60     for(ll i=0;i<m;i++)
61     {
62         ll x,y;
63         exgcd(ans,C,x,y);  
64         ll val=x*B%C;  //求丢番图方程ans*x+C*y=B的解x,即A^(m*i)*x=B(mod C)
65         val=(val%C+C)%C;  //由丢番图方程解的通式得到x的最小整数解
66         ll j=find(val);  //找到val对应的A的幂的次数
67         if(j!=-1)  return m*i+j+cnt;  //计算答案
68         ans=ans*t%C;  
69     }
70     return -1;//无解返回-1
71 }
72 void pre() {for(int i=0;i<=mod;i++)hash[i].f=0,hash[i].num=hash[i].v=-1;}
73 int main()
74 {
75     freopen("cin.in","r",stdin);
76     freopen("cout.out","w",stdout);
77     while(~scanf("%d%d%d",&A,&C,&B))
78     {
79         if(!(A+B+C))  break;
80         if(C==1)  {puts("0"); continue;}  //特判
81         pre();  A%=C;  B%=C;  //数据初始化
82         ll ans=Shank();  
83         if(ans==-1)  puts("No Solution");
84         else printf("%lld\n",ans);
85     }
86     return 0;
87 }
View Code

练习题目:【bzoj2480】Spoj3105 Mod

【bzoj3239】Discrete Logging

【bzoj2242】[SDOI2011]计算器

 

 

11、指数和原根

定义:设m>1,gcd(a,m)=1,使得 a^x mod m=1 成立的最小的 x,称为 a 对模 m 的指数(阶)。

定义:设m是正整数,a是整数,若 a 模 m 的阶等于 phi(m),则称 a 为模 m 的一个原根。

定理:如果 p 为素数,那么素数 p 一定存在原根,并且模 p 的原根的个数为 phi(p-1)。

定理:如果模 m 有原根,那么它一共有 phi ( phi (m) ) 个原根。

求模素数 p 原根的方法:对 p-1 进行素因子分解,记pi为p-1的第i个因子,若恒有a^((p-1)/pi) mod p ≠ 1  成立,则 a 就是 p 的原根。(对于合数求原根,只需把 p-1 换成 phi(p) 即可)

 1 #include<iostream>
 2 #include<cstdio>
 3 #include<cstdlib>
 4 #include<cstring>
 5 #include<ctime>
 6 #include<cmath>
 7 #include<algorithm>
 8 using namespace std;
 9 #define FILE "read"
10 int p,cnt,len,pr[50010],prime[50010],isprime[50010];
11 int gcd(int a,int b) {return !b?a:gcd(b,a%b);}
12 int fast(int a,int b,int mod) {int sum=1;for(;b;b>>=1,a=a*a%mod)if(b&1)sum*=a;return sum;}
13 inline int read()
14 {
15     int x=0,f=1;  char ch=getchar();
16     while(!isdigit(ch))  {if(ch=='-')  f=-1;  ch=getchar();}
17     while(isdigit(ch))  {x=x*10+ch-'0';  ch=getchar();}
18     return x*f;
19 }
20 void get()
21 {
22     for(int i=2;i<=50000;i++)
23     {
24         if(!isprime[i])  prime[++cnt]=i;
25         for(int j=1;j<=cnt&&prime[j]*i<=50000;j++)
26         {
27             isprime[prime[j]*i]=1;
28             if(i%prime[j]==0)  break;
29         }
30     }
31     int temp=p-1;
32     for(int i=1;i<=cnt;i++)
33     {
34         if(temp%prime[i]==0)  pr[++len]=prime[i];
35         while(temp%prime[i]==0)  temp/=prime[i];
36     }
37     if(temp>1)  pr[++len]=temp;
38 }
39 bool check(int d)
40 {
41     if(gcd(p,d)!=1)  return 0;
42     for(int i=1;i<=len;i++)  if(fast(d,(p-1)/pr[i],p)==1)  return 0;
43     return 1;
44 }
45 int main()
46 {
47     freopen(FILE".in","r",stdin);
48     freopen(FILE".out","w",stdout);
49     p=read();  get();
50     for(int i=2;i<p;i++)  if(check(i))  {printf("%d\n",i);  break;}
51     return 0;
52 }
View Code

练习题目:【POJ1284】Primitive Roots

【CF#303D】Rotatable Number 

 

 

 

参考文献——王梦迪河南省队培训讲稿《数论基础》

 

 

posted @ 2016-11-02 13:35  chty  阅读(2525)  评论(4编辑  收藏  举报