乘法逆元、扩展欧几里得算法、二元一次方程、a的n次方取余

知识点：乘法逆元，逆元的求法，二元一次方程求通解，a的n次方求余数

一，乘法逆元

---

乘法逆元的概念类似于倒数（ax=1,a−1=x），不过是在取余数的情况下的倒数。
如果(a×x)%p=1，则称x是a模p的逆元。另一种记法：ax=1( mod p)，即等式两边去膜p运算。显然x有无限多个（如果有）。

---

没有逆元我们可以很容易计算，模p的加减乘运算，但是不知道除法运算，如下所示：令,a=xp+a%p,b=yp+b%p−−−−−−−−−−−−−−−−−−−−−−−那么，a+b=(x+y)p+a%p+b%p那么，a−b=(x−y)p+a%p−b%p那么，a×b=xyp2+(x+y)p+a%p×b%p−−−−−−−−−−−−−−−−−−−−−−−所以加法，(a+b)%p=((a%p)+(b%p))%p即，a+b=(a%p)+(b%p)  ( mod p)−−−−−−−−−−−−−−−−−−−−−−−所以减法，(a−b)%p=((a%p)−(b%p))%p即，a−b=(a%p)−(b%p)  ( mod p)−−−−−−−−−−−−−−−−−−−−−−−所以乘法，(a×b)%p=((a%p)×(b%p))%p即，a×b=(a%p)×(b%p)  ( mod p)其中，x和y都是整数。

题目：http://acm.hdu.edu.cn/showproblem.php?pid=5685
逆元的作用：已知F%p和a%p的值，求(Fa)%p（我们不知道F和a的值，且F%a=0）。

---

如果已知a模p的乘法逆元为b，即(ab)%p=1那么，[Fa×(ab)]%p=[(Fa)%p×(ab)%p]%p=(Fa)%p因为，[Fa×(ab)]%p=[F×b]%p=[(F%p)×(b%p)]%p所以问题转变为求(b%p)即可

---

a模p的乘法逆元b，相当于模p运算中的a−1
模p运算中，乘以b相当于除以a，即乘以a−1。

a×b=a×a−1=1 ( mod p)[Fa]%p=[F×(1a)]%p=[F×a−1]%p=[F×b]%p[Fa×b]%p=[Fa×a−1]%p=F%p

逆元的存在性质：当a与p互素的时候，逆元才存在解，如果不互素则无解（这里说的都是整数解）。显然如果p是素数（质数），那么1到p−1的之间的数都存在模p的乘法逆元（质数和任何小于它的正整数都互质）。证明如下：

a与p不互素，则a和p存在大于1的公约数c，使得a=ca1,p=cp1。假设存在a模p的逆元b，那么，(ab)%p=1，即存在ab=xp+1，即ab−xp=1，c(a1b−xp1)=1，a1b−xp1=1c等式左边是整数，而等式右边是小数，故不存在乘法逆元。

二，求逆元

---

下面该学习如何求解乘法逆元

暴力破解：

档模p比较小的时候，我们可以枚举1到p−1，判断ax=1 ( mod p)即可。如果有逆元，那么1到p−1一定存在逆元。因为任何大于p的逆元都可以写成x=yp+x%p,由ax=ayp+a(x%p)=a(x%p)=1( modp)。所以x%p是a关于模p的逆元，且x%p小于p。

扩展欧几里得算法：

首先要了解欧几里得算法，也就是辗转相除法，gcd(a,b)=gcd(b,a%b)，每次迭代问题规模都会减小，算法复杂度是log n。扩展欧几里得算法是利用辗转相除法的逆过程来得到乘法逆元的过程。求解ax=1( mod p)，相当于求解ax=py+1，即ax+py=1，其中x和y为整数。前面说过，a和p互素才有乘法逆元，所以gcd(a,p)=1，所以相当于求解ax+py=gcd(a,p)。

---

我们辗转相除法相邻的两次迭代过程，由求gcd(a,p)转移到求gcd(p,a%p)。可以得到两个方程ax+py=gcd(a,p)=1，和px1+(a%p)y1=gcd(p,a%p)=1如果x1和y1已知的情况下，我们很容易得到一组x和y的解，因为，ax+py=gcd(a,p)=1=gcd(p,a%p)=px1+(a%p)y1，得到ax+py=px1+(a%p)y1=px1+(a−(a/p)∗p)y1=ay1+p(x1−(a/p)y1)即，x=y1，y=x1−(a/p)y1。所以求解方程ax+py=gcd(a,p)的解，转化为求解方程px1+(a%p)y1=gcd(p,a%p)的解。和辗转相除法的过程是一样的，每次迭代，问题的规模都在减小。最后会化简为求解gcd(a,p)xn+0∗yn=gcd(a,p)，。这个时候只需要令xn=1，yn=0即可，倒退回去可以求得x和y的值。

---

上述过程求解ax+py=gcd(a,p)的过程，并不要求gcd(a,p)一定等于1。由上述过程 ，我们可以推广得到求解任意二元一次方程ax+by=c的解。

显然可以先求解ax′+by′0=gcd(a,b)，然后得到原方程的一组解x0=x′∗(cgcd(a,b))，y0=y′∗(cgcd(a,b))。然后通解的表达式为x=x0+bgcd(a,b)∗t，y=y0−agcd(a,b)∗t，其中t为任意整数。因为显然ax+by=ax0+by0+a∗bgcd(a,b)∗t−b∗agcd(a,b)∗t=ax0+by0。

需要注意的一点是只有当c%gcd(a,b)=0时，方程ax+by=c才有整数解证明如下：

因为显然可以ax+by=gcd(a,b)∗(a′x+b′y)=c，如果x和y都是整数，而a′和b′也是整数，所以c必须是gcd(a,b)的倍数，才有整数解。

---

扩展欧几里得算法代码如下：

int exgcd(int a, int b, int &x, int &y){    if(b==0){        x = 1;        y = 0;        return a;    }else{        int _gcd = exgcd(b,b%a,x,y);        int tmp = x - (a/b)*y;        x = y;        y = tmp;        return _gcd;    }}x = (x % p + p) % p; // 通常需要求解最小的正整数逆元解ax=1(mod p)// 再次提醒a和p互素，即gcd(a,p) = 1

费马小定理：

费马小定理：假如p是质数，且gcd(a,p)=1，那么ap−1≡1 ( mod p)，即：假如a是整数，p是质数，且a,p互质(即两者只有一个公约数1)，那么a的(p-1)次方除以p的余数恒等于1。

由ap−1=1 ( mod p)很容易得到a∗ap−2=1 ( mod p)，那也就是a的逆元是ap−2。

需要注意的是需要判断a和p是否互素，如果用gcd来判断的话，那还不如直接用扩展欧几里得算法。

O(n)求1~n逆元表：

由前面可知，当p为素数时，1到p−1都存在模p的乘法逆元。
采用递推的方法来求解所有1到p−1的逆元。如下：

显然当a=1的时候，a的逆元是1。当a大于1的时候，由p%a+a∗(p/a)=p，得到−a∗(p/a)=p%a ( mod p)，等式两边都乘以p%a的逆元(p%a)−1，得到−a∗(p/a)∗(p%a)−1=1 ( mod p)，即：a∗[−(p/a)∗(p%a)−1]=1 ( mod p)，所以a的逆元是−(p/a)∗(p%a)−1，因为p%a小于a，所以(p%a)−1是已知的，所以从1到p−1的顺序依次计算出所有数字的逆元。

也注意上面求出来的逆元是负数，很多时候需要转化为正整数，而且是最小的正整数逆元。

为了得到正数，由a∗[−(p/a)∗(p%a)−1]=1 ( mod p)，得到a∗[[p−(p/a)]∗(p%a)−1]=1 ( mod p)，再由[[p−(p/a)]∗(p%a)−1]%p得到最小的正整数。

代码如下：

const int p = 13;int inv[p+2];inv[1] = 1;for (int i = 2;i < p;i++)    inv[i] = ((p - p / i)*inv[p%i]) % p;

对于1到p的都只存在一个小于p的逆元，为什么？

首先a和p的最大公约数gcd(a,p)是1，求逆元相当于，求解方程ax+py=1，由前面的通解表达式可得到x=x0+pgcd(a,p)t=x0+pt，可以看到，通解x就是p的等差数列，所以只能存在一个。

三，a的n次方取余

---

如求20172017%1777或者2100%9973的值：

前面我们知道a×b=(a%p)×(b%p)  ( mod p)所以这里二分法分解原式子，将问题规模变小：20172017%1777={20171008×20171008×2017}%1777，显然问题化简为了求20171008%1777，进而可以化简为20171008%1777={2017504×2017504×1}%1777容易总结出规律，问题规模在成指数级别减小。递归求解公式为：ab%p={(ab/2%p)2×ab%2}%p，时间复杂度是log n。

代码如下：

long long remainder_of_an_exponential_recur(int a, int b, int c){  // 递归求 (a^b)%c    if (b == 0)        return 1;    long long tmp = remainder_of_an_exponential_recur(a, b >> 1, c);    if (b & 1)        return (tmp*tmp*a) % c;    else        return (tmp*tmp) % c;}long long remainder_of_an_exponential_loop(int a, int b, int c){  // 递推求 (a^b)%c    long long res = 1;    long long tmp = a;    while (b)    {        if (b & 1) res = (res*tmp) % c;        tmp = (tmp*tmp) % c;        b >>= 1;    }    return res;}

---

如果除数是2i，求余数。这个比较简单a%2也就是求a的二进制表示方式中的最后一位的值。a%2=a%(21)=a&20=a&1。a%8也就是求a的二进制表示方式中的最后三位的值。a%8=a%(23)=a&(23−1)=a&7。以此类推a%(2i)=a&(2i−1)。

---

最后一个需要说的是公式：ab%p=a%(b×p)b%p。（其中a%b=0）。这个公式一般都是在无法计算得到a的情况下使用的，这个公式当然也可以两边乘以b的逆元。这个公式画图理解比较好，如下所示：

在ALL X问题中：http://acm.hdu.edu.cn/showproblem.php?pid=5690。
要求检验一个全由x组成的m位数字F(x,m)是否可以满足等式F(x,m) mod k=c。

既可以利用余数有限重复循环的方法，即按照除法的方法依次计算余数，如果余数重复，那么就开始进入了循环，而余数的个数是有限的，小于k个。这和求a/b的小数表示形式是一样的（它是无限的添加0，本题中是无限的添加x）。

另一种方法是求解F(x,m)，容易知道F(x,m)=x×10m−19，那么F(x,m) mod k=[(x%k)×10m−19%k]%k，相当于求解10m−19%k，利用上述公式可以得到10m−19%k=(10m−1)%(9k)9%k=[10m%(9k)]−19%k。而[10m%(9k)]可以利用a的n次方取余的方法来求解。

求an，可以用快速幂乘，a可以是矩阵。
an=an/2∗an/2,n是偶数
a^{n}=a^{n/2}*a^{n/2}*ahttps://www.baidu.com/s?ie=UTF-8&wd=快速幂乘&tn=98012088_4_dg&ch=10  
,n是奇数数<script type="math/tex" id="MathJax-Element-1408">a^{n}=a^{n/2}*a^{n/2}*ahttps://www.baidu.com/s?ie=UTF-8&wd=快速幂乘&tn=98012088_4_dg&ch=10 ,n是奇数数</script>

参考：
http://blog.csdn.NET/tsaid/article/details/7365936
http://www.cnblogs.com/ka200812/archive/2011/09/02/2164404.html
https://wenku.baidu.com/view/c1b06ea60029bd64783e2c63.html
http://blog.csdn.net/stcyclone/article/details/52081822