1、快速幂
计算a^b的快速算法，例如，3^5，我们把5写成二进制101，3^5=3^1*1+3^2*2+3^4*1
ll fast(ll a,ll b){ll ans=;for(;b;b>>=,a=mul(a,a))if(b&)ans=mul(ans,a);return ans;}//一行快速幂
2、快速乘
当模数较大时，直接乘会爆掉long long,需要快速乘法。
即用浮点计算倍数，做差相当于计算余数模2^63的结果，然后再模一下就好了（因为余数不超过long long）
typedef long long ll;
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)的公约数是一样的，其最大公约数也必然相等，得证。
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
扩展欧几里得的过程可以理解为从末状态向上不断回溯的过程，直到得到原方程的一组解。
void exgcd(int a,int b,int &x,int &y)
{
if(b==) {x=; y=; return;}
exgcd(b,a%b,x,y);
int t=x;x=y;y=t-a/b*y;
}
丢番图方程的通解：若(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的乘法逆元。）
void China()
{
for(ll i=;i<=n;i++)
{
ll x,y,Mi=M/m[i];
exgcd(Mi,m[i],x,y);
ans=(ans+Mi*x*a[i])%M;
}
}
<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
int China()
{
int A=a[],k,y; M=m[];
for(int i=;i<=n;i++)
{
int da=a[i]-A,g;
exgcd(M,m[i],g,k,y);
if(da%g) return -;
k*=da/g;
A+=k*M;
M=M*m[i]/g;
A=(A+M)%M;
}
return A;
}
练习题目：【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)
void get()
{
memset(check,,sizeof(check));
int cnt=;
for(int i=;i<=N;i++)
{
if(!check[i]) {prime[++cnt]=i; phi[i]=i-;}
for(int j=;j<=cnt&&prime[j]*i<=N;j++)
{
check[i*prime[j]]=;
if(i%prime[j]) phi[i*prime[j]]=phi[i]*(prime[j]-);
else {phi[i*prime[j]]=phi[i]*prime[j]; break;}
}
}
}
练习题目：【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的所有素数。
int cnt,prime[MAXN],check[MAXN];
void get()
{
for(int i=;i<=n;i++)
{
if(!check[i]) prime[++cnt]=i;
for(int j=;j<=cnt&&prime[j]*i<=n;j++)
{
check[prime[j]*i]=;
if(i%prime[j]==) break;//保证线筛的时间复杂度
}
}
}
<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的个数
int fast(int a,int b,int mod) {int sum=;for(;b;b>>=,a=a*a%mod)if(b&)sum=sum*a%mod;return sum;}//一行快速幂
bool Rabin_Miller(int p,int a)
{
if(p==) return ;
if(p&==||p==) return ;
int d=p-;
while(!(d&)) d>>=;
int m=fast(a,d,p);
if(m==) return ;
for(;d<p;d<<=,m=m*m%p) if(m==p-) return ;
return ;
}
bool isprime(int x)
{
static int prime[]={,,,,,,,,};//相当于全局变量
for(int i=;i<;i++)
{
if(x==prime[i]) return ;
if(!Rabin_Miller(x,prime[i])) return ;
}
return ;
}
练习题目：【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))
int cnt,fac[],num[];
int gcd(int a,int b) {return !b?a:gcd(b,a%b);}
int fast(int a,int b,int mod) {int sum=;for(;b;b>>=,a=a*a%mod)if(b&)sum=sum*a%mod;return sum;}//一行快速幂
bool Rabin_Miller(int p,int a)
{
if(p==) return ;
if(p&==||p==) return ;
int d=p-;
while(!(d&)) d>>=;
int m=fast(a,d,p);
if(m==) return ;
for(;d<p;d<<=,m=m*m%p) if(m==p-) return ;
return ;
}
bool isprime(int x)
{
static int prime[]={,,,,,,,,};//相当于全局变量
for(int i=;i<;i++)
{
if(x==prime[i]) return ;
if(!Rabin_Miller(x,prime[i])) return ;
}
return ;
}
void Pollord_Rho(int x)
{
if(isprime(x)) {fac[++cnt]=x; return;}
int c=;
while()
{
int x1(),x2(),k(),i();
while()
{
x1=((x1*x1)%x+c)%x;
int d=gcd(abs(x1-x2),x);
if(d>&&d<x)
{
Pollord_Rho(d);
Pollord_Rho(x/d);
return;
}
if(x1==x2) break;
if(++i==k) k<<=,x2=x1;
}
c++;
}
}
void solve(int x)
{
Pollord_Rho(x);
sort(fac+,fac+cnt+);
int len=;
for(int i=;i<=cnt;i++)
{
if(fac[i]==fac[i-]) num[len]++;
else num[++len]=,fac[len]=fac[i];
}
for(int i=;i<len;i++) printf("%d^%d*",fac[i],num[i]);
printf("%d^%d\n",fac[len],num[len]);
}
练习题目：【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的代码作为模板（良心的博主加上了注释）：
/*************
bzoj 2480
by chty
2016.11.8
*************/
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<ctime>
#include<cmath>
#include<algorithm>
using namespace std;
#define mod 99991
typedef long long ll;
struct node{ll v,num,f;}hash[mod+];
ll A,B,C;
inline ll read()
{
ll x=,f=; char ch=getchar();
while(!isdigit(ch)) {if(ch=='-') f=-; ch=getchar();}
while(isdigit(ch)) {x=x*+ch-''; ch=getchar();}
return x*f;
}
ll gcd(ll a,ll b) {return !b?a:gcd(b,a%b);}
void insert(ll v,ll x) //哈希表插入操作
{
ll t=v%mod;
while(hash[t].f&&hash[t].v!=v) {t++; if(t>mod) t-=mod;}
if(!hash[t].f) {hash[t].f=;hash[t].num=x;hash[t].v=v;}
}
ll find(ll v) //哈希表查询操作
{
ll t=v%mod;
while(hash[t].f&&hash[t].v!=v) {t++; if(t>mod) t-=mod;}
if(!hash[t].f) return -;
else return hash[t].num;
}
void exgcd(ll a,ll b,ll &x,ll &y) //扩展欧几里得算法求丢番图方程的一组解
{
if(!b) {x=; y=; return;}
exgcd(b,a%b,x,y);
ll t=x;x=y;y=t-a/b*y;
}
ll Shank()
{
ll ret=;
for(ll i=;i<=;i++) {if(ret==B)return i;ret=ret*A%C;}
//判断在0~50内是否有解，这就是Baby_Step，不加这一步会wa掉，因为下面做除法时，假设的是A足够用
ll temp,ans(),cnt();
while((temp=gcd(A,C))!=)//当A与C不互质时，拿出若干个A出来与C来消去公共因子，使得gcd(A,C)=1为止
{
if(B%temp) return -;//根据裴蜀定理，如果gcd(A,C)|B不成立，则无解
C/=temp; B/=temp;
ans=ans*(A/temp)%C;
cnt++; //记录已经乘了cnt个A，以及这些A与C消掉公因子后的结果ans
}
ll m=(ll)ceil(sqrt(C*1.0)),t(); //保证时间复杂度为O(n^0.5)
for(ll i=;i<m;i++) {insert(t,i);t=t*A%C;} //把A^0~A^m存到哈希表中
for(ll i=;i<m;i++)
{
ll x,y;
exgcd(ans,C,x,y);
ll val=x*B%C; //求丢番图方程ans*x+C*y=B的解x，即A^(m*i)*x=B(mod C)
val=(val%C+C)%C; //由丢番图方程解的通式得到x的最小整数解
ll j=find(val); //找到val对应的A的幂的次数
if(j!=-) return m*i+j+cnt; //计算答案
ans=ans*t%C;
}
return -;//无解返回-1
}
void pre() {for(int i=;i<=mod;i++)hash[i].f=,hash[i].num=hash[i].v=-;}
int main()
{
freopen("cin.in","r",stdin);
freopen("cout.out","w",stdout);
while(~scanf("%d%d%d",&A,&C,&B))
{
if(!(A+B+C)) break;
if(C==) {puts(""); continue;} //特判
pre(); A%=C; B%=C; //数据初始化
ll ans=Shank();
if(ans==-) puts("No Solution");
else printf("%lld\n",ans);
}
return ;
}
练习题目：【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) 即可）
#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<ctime>
#include<cmath>
#include<algorithm>
using namespace std;
#define FILE "read"
int p,cnt,len,pr[],prime[],isprime[];
int gcd(int a,int b) {return !b?a:gcd(b,a%b);}
int fast(int a,int b,int mod) {int sum=;for(;b;b>>=,a=a*a%mod)if(b&)sum*=a;return sum;}
inline int read()
{
int x=,f=; char ch=getchar();
while(!isdigit(ch)) {if(ch=='-') f=-; ch=getchar();}
while(isdigit(ch)) {x=x*+ch-''; ch=getchar();}
return x*f;
}
void get()
{
for(int i=;i<=;i++)
{
if(!isprime[i]) prime[++cnt]=i;
for(int j=;j<=cnt&&prime[j]*i<=;j++)
{
isprime[prime[j]*i]=;
if(i%prime[j]==) break;
}
}
int temp=p-;
for(int i=;i<=cnt;i++)
{
if(temp%prime[i]==) pr[++len]=prime[i];
while(temp%prime[i]==) temp/=prime[i];
}
if(temp>) pr[++len]=temp;
}
bool check(int d)
{
if(gcd(p,d)!=) return ;
for(int i=;i<=len;i++) if(fast(d,(p-)/pr[i],p)==) return ;
return ;
}
int main()
{
freopen(FILE".in","r",stdin);
freopen(FILE".out","w",stdout);
p=read(); get();
for(int i=;i<p;i++) if(check(i)) {printf("%d\n",i); break;}
return ;
}
练习题目：【POJ1284】Primitive Roots
【CF#303D】Rotatable Number
参考文献——王梦迪河南省队培训讲稿《数论基础》