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;}//一行快速冪
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;}//一行快速乘
練習題目：【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
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 }
丟番圖方程的通解：若(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 }
<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 }
練習題目：【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 }
練習題目：【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 }
<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 }
練習題目：【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 }
練習題目：【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 }
練習題目：【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 }
練習題目：【POJ1284】Primitive Roots
【CF#303D】Rotatable Number
參考文獻——王夢迪河南省隊培訓講稿《數論基礎》