常用的数论算法

最大公约数与最小公倍数

int gcd(int a, int b)
{
    b==0?a:gcd(b,a%b);
}

int lcm(int a,int b)
{
    return a/gcd(a,b)*b;
}

---

二进制运算

1.按位异或^

功能：	当两个数的二进制表示，进行异或运算时，当前位的两个二进制表示不同则为1相同则为0.同为0异为1
特点：
	(1) 0^0=0,0^1=1 0异或任何数＝任何数
	(2) 任何数异或自己＝把自己置0

2.位运算

与（&）	进制位全为1才为1	0 & 0 = 0	1 & 0 = 0	0 & 1 = 0	1 & 1 = 1
或（|）	进制位全为0才为0	0 | 0 = 0	1 | 0 = 1	0| 1 = 1	1 | 1 = 1
异或（^）	同为0异为1	0 ^ 0 = 0	1 ^ 0 = 1	0 ^ 1 = 1	1 ^ 1 = 0

记忆方法:	使二进制中的1与逻辑运算中的true对应，0与逻辑运算中的false对应
	由真值表可知：对于&运算全为真才为真，对于|运算有真则为真全为假才为假。

用途：	&运算判断奇偶：a&1==1则a是奇数,a&1==0则a是偶数
	|运算实现偶数的自增1：a是偶数则a|1==a+1;a是奇数a|1==a;任何数|0==自身

3.位移运算

左移运算：

左移运算符m<<n表示吧m左移n位。左移n位的时候，最左边的n位将被丢弃，同时在最右边补上n个0.比如：
00001010 << 2 = 00101000

10001010 << 3 = 01010000
常用的是实现乘法运算a<<1==a*2，a<<2=a*4;

右移运算

右移运算符m>>n表示把m右移n位。右移n位的时候，最右边的n位将被丢弃。但右移时处理最左边位的情形要稍微复杂一点。这里要特别注意，如果数字是一个无符号数值，则用0填补最左边的n位。如果数字是一个有符号数值，则用数字的符号位填补最左边的n位。也就是说如果数字原先是一个正数，则右移之后再最左边补n个0；如果数字原先是负数，则右移之后在最左边补n个1.下面是堆两个8位有符号数作右移的例子：

00001010 >> 2 = 00000010

10001010 >> 3 = 11110001

---

快速幂O(log₂N)

思路：假设我们要求a^b，那么其实b是可以拆成二进制的，该二进制数第i位的权为2^(i-1)，例如当b==11时a^11=a^(2^0+2^1+2^3)11的二进制是1011，11 = 2³×1 + 2²×0 + 2¹×1 + 2º×1，因此，我们将a¹¹转化为算 a^(2^0)a^(2^1)*a^(2^3)=a^1*a^2 a^8，看出来快的多了吧原来算11次，现在算三次。

注意a的幂每次都是2^n(n=0,1,2,3……)递增1，2，4，8,…….

//根据实际情况改变返回类型，以及是否对某个数取余
int qpow(int a,int b)
{
    int ans = 1,base = a;
    while(b!=0)
    {
        if(b&1)//二进制末位为1就乘起来
            ans *= base;
        base *= base;//a的幂的自增
        b>>=1;
    }
    return ans;
}

程序过程： b=11（1011）

初始:	ans=1	base=a	b=1011
b&1=1	ans=a	base=a^2	b=101
b&1=1	ans=a*a^2	base=a^4	b=10
b&1=0	ans=a*a^2	base=a^8	base=1
b&1=1	ans=a*a^2*a^8	base=a^16	base=0

结果就是：ans=a*a^2*a^8

调试快速幂的方法：随便找1个质数p，检测a^(p-1)%p是否都为1就行啦 (1<=a＜p)

---

素数快速判断

bool is_Prime( int num )
{
    //快速判断一个数是否为素数

    //两个较小数另外处理
    if(num ==2|| num==3 )
        return 1 ;
    //不在6的倍数两侧的一定不是质数
    if(num %6!= 1&&num %6!= 5)
        return 0 ;
    int tmp =sqrt( num);
    //在6的倍数两侧的也可能不是质数
    for(int i= 5; i <=tmp; i+=6 )
        if(num %i== 0||num %(i+ 2)==0 )
            return 0 ;
    //排除所有，剩余的是质数
    return 1 ;
}

埃氏筛素数打表

找出n以内的所有素数：做法：做法其实很简单，首先将2到n范围内的整数认定为素数，其中2是最小的素数。将表中所有的2的倍数划去，表中剩下的最小的数字就是3，他不能被更小的数整除，所以3是素数。再将表中所有的3的倍数划去……以此类推，如果表中剩余的最小的数是m，那么m就是素数。然后将表中所有m的倍数划去，像这样反复操作，就能依次枚举n以内的素数，这样的时间复杂度是O(nloglogn)。

const int MAXN=1000;
int prime[MAXN];//记录素数
bool is_pri[MAXN + 10]; //is_pri[i]表示i是否为素数
//返回n以内素数的个数
int sieve(int n)
{
    int p = 0;
    for(int i = 0; i <= n; i++)
        is_pri[i] = true; //初始化全部是素数
    is_pri[0] = is_pri[1] = false;//0和1不是素数
    for(int i = 2; i <= n; i++)
    {
        if(is_pri[i])
        {
            prime[++p] = i;
            for(int j = 2 * i; j <= n; j += i)
                if(is_pri[j])
                is_pri[j] = false;//i的倍数全部剔除
        }
    }
    return p;
}

扩展欧几里德算法

贝祖（裴蜀）等式

数论中，贝祖定理是一个关于最大公约数gcd（）的定理，说明了对于任何正整数a、b和他们的最大公约数d，关于未知数x和y的线性丢番图方程（称为贝祖等式）：
$$ax+by=m有整数解当且仅当m是d的倍数$$
贝祖等式有解时必然有无穷组解（x，y），可用扩展欧几里德算法求的一组特解
特别的 $$ax+by=1有解当且仅当a，b互质$$

---

求解ax+by=gcd（a,b）

扩展欧几里得算法就是在求 a,b 的最大公约数 m=gcd(a,b) 的同时，求出贝祖等式ax + by = m的一个特解 (x,y)。

typedef long long ll;
void exgcd(ll a, ll b, ll& d, ll& x, ll &y)
{//d是gcd（a,b）  （x，y）是一组特解
    if(!b)
    {
        d = a, x = 1, y = 0;
    }
    else
    {
        exgcd(b, a % b, d, y, x);
        y -= x * (a / b);
    }
}

求解不定方程

$$求解ax+by=c的通解（等式有解当且仅当c\%gcd（a，b）==0）。如果有解，则令 K = c/d；$$
$$那么，对方程aX+bX = d，假设由拓展欧几里得求出一组解为(X0,Y0)，那么有aX0+bY0 = d；$$
$$等式两边同时乘以K，即K*( aX0+bY0 ) = d*K = c；$$

$$由恒等关系，原方程的解（x0，y0）： x0 = KX0 = c/d * X0，y0 = KY0 = c/d *Y0。$$
$$不定方程的通解： 若（x0，y0）是不定整数方程ax+by = c的一组解，则他的任意整数解都可以表示成（x0+ tb’, y0-ta’）t为常数.$$ $$其中a’ = a/gcd(a,b), b’ = b/gcd(a,b).$$

exgcd(a,b,d,X0,Y0);//d,x,y均为引用传参 
if(c%d==0)//满足贝祖等式的条件
   K=c/d;
 else
    cout<<"无解";
  int x=KX0,y=KY0;//不定方程一组特解
  int A=a/d,B=b/d
  for(int i=1;;i++)
  {
     x+=B;//通解
     y+=A;
  }

欧拉函数

欧拉函数phi(x)是指不大于正整数x的与x互质的正整数的个数。例如phi(1)=1，phi(2)=1，phi(3)=2，phi(4)=2，phi(5)=4，phi(6)=2等等。很显然，对每一个质数p，phi(p)=p-1。而对每一个质数的幂phi(p^n)=(p-1)×p^(n-1)。欧拉函数是积性函数，如果m、n互质，那么phi(m×n)=phi(m)×phi(n)。
求欧拉函数的值需要用到算术基本定理。将n写成其质因子幂的连乘积的形式，n=p1^k1×p2^k2×…×pr^kr。那么，phi(n)=n×(1-1/p1)×…×(1-1/pr)。利用这个公式可以写出求欧拉函数的代码。当然，首先要求质数（素数筛选法）。如果问题范围是MAXSIZE，那么只需求出SIZE=sqrt(MAXSIZE)以内的质数即可。即保证任何一个数的最小质因子一定在这其中。

//先用素数筛选法求出sqrt（MAXSIZE）以内的质数
int euler(int n){
    int ret = n;
    for(int i=0;P[i]*P[i]<=n;++i){ //P是保存质数的数组
        if ( n % P[i] ) continue;

        ret -= ret / P[i];//计算欧拉函数
        n /= P[i];
        while( 0 == n % P[i] ) n /= P[i];
        if ( 1 == n ) break;//这三行用于加速并有助于最后得到大的质因子
    }
    if ( n != 1 ) ret -= ret / n;  //假如n含有超过SIZE的质因子
    return ret;
}

逆元（inv）

什么是逆元

定义:对于正整数a和m，使得

$$a⋅x≡1(mod \ \ m)$$ 成立的x的 最小正整数解，是a关于m的逆元。
$$欧拉定理：a^{phi(m)}=1\ \ (mod\ \ m) \ $$
$$由定义：a⋅x≡1\ \ (mod \ m)$$
$$可知x=a^{phi(m)-1}\ \ \ \ \ \ 也可知有a^{-1}=a^{phi(m)-1}由此可知逆元使得除法变为乘法$$

逆元的作用

当求解公式：(a/b) （mod）m 时，因b可能会过大，会出现爆精度的情况，所以需变除法为乘法：

$$\ 设c是b关于m的逆元，则有b*c≡1(mod m)\ \ \ \ \ \ \ \ \ 即c=b^{-1}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ \ \ \ \ \ \ \ \ $$

$$\frac{a}{b}\ (mod\ \ m) = \frac{a}{b}\ *1(mod\ \ m) = \frac{a}{b}\ *b*c(mod\ \ m) = a*c(mod \ \ m);$$

$$结论：为求\frac{a}{b}\ (mod\ \ m)（b较大会爆精度），已知b的逆元是c，则有 \frac{a}{b}\ (mod\ \ m)=a*c(mod\ \ m)。$$ 那么求解 问题的关键就是求解b（mod）m的逆元c。

---

逆元的求法

利用欧拉函数对模m是否为质数没要求

$$由欧拉定理可知:x=a^{phi(m)-1}是a关于m的逆元\ \ \ \ \ 且phi（m）可用欧拉函数求出$$

x=euler(m);//欧拉函数求出phi（m）
ans=qpow(a,x-1);//快速幂求a^(phi(m)-1)

费马小定理

在 模m 是素数的情况下，对任意整数 a 都有 $$a^m≡a(mod\ \ m)$$

如果 x 无法被m 整除（x与m互质），则有

$$a^{m−1} ≡ 1(mod\ \ m)$$
可以在 m 为素数的情况下求出一个数的逆元， $$a∗a^{m-2}≡1(mod\ \ m) \ \ \ \ \ \ 则a^{m−2}即为a关于m的逆元。$$
$$由欧拉函数性质:对每一个质数m，phi(m)=m-1\ \ \ \ \ 以及欧拉定理:a^{-1}=x^{phi(m)-1}可知费马小定理是欧拉定理的特例$$

//用快速幂实现x^(m-2)
const int mod = 1e9+9;//可以根据a的数据范围而定mod>a且为素数由质数的性质:一个质数，另一个不为它的倍数，这两个数为互质数,可证逆元存在
long long qpow(long a,long b)
｛
    if(b<0) return 0;
    long long ans=1,base=a;
    while(b)
    {
       if(b&1)
           ans=(ans*base)%mod;//放止溢出
        base=(ans*base)%mod;
        b>>=1;
     }
     return ans;
 ｝
long long inv(long long a)
{
    return qpow(a,mod-2);
 }
 int check(long long ans,long long a)
 {//检验a是否有逆元
     if(ans*a)==1
       return 1;
     return 0;
  }

扩展欧几里德求解逆元

给定模数m，求a的逆相当于求解ax=1(mod m)
这个方程可以转化为ax-my=1
然后套用求二元一次方程的方法，用扩展欧几里得算法求得一组x0,y0和gcd
检查gcd是否为1
gcd不为1则说明逆元不存在 （逆元存在条件a与m互质）
若为1，则调整x0到0~m-1的范围中即可
PS：这种算法效率较高，常数较小，时间复杂度为O(ln n)

typedef  long long ll;
void extgcd(ll a,ll b,ll& d,ll& x,ll& y){
    if(!b)
    {
      d=a;
      x=1; 
      y=0;
      }
    else
    {
     extgcd(b,a%b,d,y,x); 
     y-=x*(a/b); 
    }
}
ll inverse(ll a,ll n){
    ll d,x,y;
    extgcd(a,n,d,x,y);
    return d==1?(x+n)%n:-1;
}

逆元的调试：若ans*a==1答案正确

逆元打表与数据规模线性相关

有时会遇到这样一种问题，在模质数p下，求1~n逆元 n< p（这里为奇质数）。可以O(n)求出所有逆元，有一个递推式如下

$$inv[i]=(m-m/i)*inv[M\%i]\%M$$

typedef  long long ll;
const int N = 1e5 + 5;
int inv[N];

void inverse(int n, int p) {
    inv[1] = 1;
    for (int i=2; i<=n; ++i) {
        inv[i] = (ll) (p - p / i) * inv[p%i] % p;
    }
}

上面说的那么多是为了求解b（mod m）的逆元，以便求解 $$\frac{a}{b}\ (mod\ \ m)$$ 然而求解逆元的前提b与m互质。问题来了如果b与m不互质怎么求 $$\frac{a}{b}\ (mod\ \ m)$$

$$ans=\frac{a}{b}\ (mod\ \ m)=\frac{a(mod(m*b))}{b}$$

PS:实际上$\frac{a}{b} mod (b*m)$这种的对于所有的都适用，不区分b与m互不互素,而费马小定理和扩展欧几里得算法求逆元是有局限性的，它们都会要求与互素，如果a与m不互素，那就没有逆元，这个时候需要$\frac{a}{b} mod (b*m)$来搞（此时就不是逆元的概念了）。但是当a与m互素的时候，bm可能会很大，不适合套这个一般公式，所以大部分时候还是用逆元来搞

===

总结

逆元的知识有点多梳理一下

$$首先逆元是为了解决\frac{a}{b}\ (mod\ \ m)的搞精度问题$$

逆元的存在条件b与m互质

求解方法：	适用范围	优缺点
欧拉定理 $$b^{phi(m)-1}是b关于m的逆元$$	只要b与m互质即可	涉及到素数打表、快速幂、欧拉函数求解问题（不常用）
扩展欧几里德	只要b与m互质	只涉及到扩展欧几里德算法（方法比欧拉定理简单，适用范围比费马大）这种算法效率较高，常数较小，时间复杂度为O(ln n)
费马小定理 $$b^{m-2}即为所求$$	要求b与m互质且m是质数ps:若m是质数且a不是m的倍数则a与m互质	只涉及到快速幂的求解这种算法复杂度为O(log2 n) 在几次测试中，常数似乎较上种方法大
逆元打表 $$inv[i]=(m-m/i)*inv[M\%i]\%M$$	在模质数p下，求1~n逆元 n< p（这里为奇质数）打表方法	复杂度O（n）

若b与m不互质不能用逆元去求用通式： $$ans=\frac{a}{b}\ (mod\ \ m)=\frac{a(mod(mb))}{b}$$
缺点是m*b如果太大会溢出

---

求组合数

打表：适用范围maxn=1003

$$C(n,k)= C(n−1,k−1)+C(n−1,k)$$
$$C(n,0)=C(n,n)= 1$$

const int maxn=1005;
const int mod=100003;
int ZHS[maxn][maxn];

void calcC(int n)
{ //打表n<1005
    for(int i=0;i<n;i++)
    {
        ZHS[i][0]=1;
        for(int j=1;j<i;j++)
            ZHS[i][j]=(ZHS[i-1][j-1]+ZHS[i-1][j])%mod;
        ZHS[i][i]=1;
    }
}

int main()
{
  calcC(1005);
  int a,b;
  while(cin>>a>>b)
  cout<<ZHS[a][b]<<endl;
  return 0;
}

乘法逆元求组合数取模

$$较大的组合数一般会对素数取模：C(n, r) = \frac {n!} {r! * (n - r)!}$$ $$而我们的目标就是计算 C(n,r)\%\ p =n！*inv（r！*（n-r）！））\ \%p的值$$ $$\ \ p一般为质数 且1\leq n、m、p\ \leq10^5由之前对逆元的知识知道我们可以利用逆元求解。$$

//求逆元的方法很多（详情看上面关于逆元的知识）
//由于p一般为质数...所以费马小定理、扩展欧几里德、欧拉公式、线性打表
//甚至求模的通式都可以使用（注意是否会爆数据就好）
int inv(int a) 
{//线性求逆元公式简单   
    return a == 1 ? 1 : (long long)(MOD - MOD / a) * inv(MOD % a) % MOD;  
} 
LL comb(LL n,LL m)  
{  
    if(m < 0)return 0;  
    if(n < m)return 0;  
    if(m > n-m) m = n-m;  

    LL up = 1, down = 1;  
    for(LL i = 0 ; i < m ; i ++){  
        up = up * (n-i) % MOD;//求分子与分母的阶乘  
        down = down * (i+1) % MOD;  
    }  
    return up * inv(down) % MOD;  
}  

Lucas定理求大组合数求模

$$定义:Lucas(n,m,p)=C(n,m)\ \%p$$
$$Lucas定理：Lucas(n,m,p)=C(n\%p,m\%p)*Lucas(n/p,m/p,p)$$
$$Lucas定理使得组合数的均参数<p，而且p\leq10^5。那么就可以通过逆元求解组合数，在通过Lucas定理优化参数进而求解出大的组合数模p的值$$

 LL Lucas(LL n, LL m, int p)
 {
         return m ? Lucas(n/p, m/p, p) * comb(n%p, m%p, p) % p : 1;       
 }

总结

组合数求解方法：	范围
暴力求解 $$C(n, r) = \frac {n!} {r! * (n - r)!}$$	很小
公式打表 $$C(n,k)= C(n−1,k−1)+C(n−1,k) \ \ \ 、C(n,0)=C(n,n)= 1$$	$$10^3左右$$
逆元求解 $$C(n,r)\%\ p =n！*inv（r！*（n-r）！））\ \%p$$	$$n、m\leq10^9、p\leq10^5$$
$$Lucas定理：Lucas(n,m,p)=C(n\%p,m\%p)*Lucas(n/p,m/p,p)$$	$$p\leq10^5,n、m\leq10^{18}$$