第十一课：数学基础（快速幂、逆元、扩展欧几里得算法、中国剩余定理）

BlessingSoftware1

已于 2022-09-22 08:38:25 修改

阅读量634

点赞数

分类专栏：摆烂日记from acwing 文章标签：算法

于 2022-07-08 23:11:04 首次发布

本文链接：https://blog.csdn.net/weixin_64197656/article/details/125667128

版权

摆烂日记from acwing 专栏收录该内容

45 篇文章 1 订阅

订阅专栏

〇、欧拉定理

若a与N互质，则a的phi（N）次方模N同余1。

若N是一个质数，则phi（N）=N-1：（费马小定理）

一、快速幂

（1）定义及求法

用途：快速求出a的k次方模p的值。

方式：1.预处理出a的2的x次方模p的值的序列。其中x最大到log（k）。

则:

$gif.latex?a%5E%7Bk%7D%3Da%5E%7B2%5E%7Bx_%7B1%7D%7D%7D*a%5E%7B2%5E%7Bx_%7B%7D%7D%7D...*a%5E%7B2%5E%7Bx_%7Bn%7D%7D%7D$

$gif.latex?k%3D2%5E%7Bx_%7B1%7D%7D+2%5E%7Bx_%7B2%7D%7D+...+2%5E%7Bx_%7Bn%7D%7D$

此时，发现k用二进制表示出来则可相应求得x的所有值。因此可以替换：

$gif.latex?a%5E%7Bk%7Dmod%28p%29%3Da%5E%7B2%5E%7Bx_%7B1%7D%7D%7D*a%5E%7B2%5E%7Bx_%7B%7D%7D%7D...*a%5E%7B2%5E%7Bx_%7Bn%7D%7D%7D%20mod%28p%29$

至于预处理这个序列，发现a的2的0次方a，a的2的一次方为a的平方。实际上后一项就是前一项的平方。因此预处理也比较方便。

eg：

快速幂的实际代码实现起来，融合了预处理与计算，需要牢记：

//这里填你的代码^^
#include<iostream>
#include<algorithm>
using namespace std;

typedef long long LL;


LL qmi(int a,int b,int p)
{
    //考虑b的二进制形式
    LL res=1;
    while(b)
    {
        if(b&1) res=res*a%p;//如果b该位有1
        b>>=1;//去除一位
        a=a*(LL)a%p;//a平方
    }
    return res;
}


int main()
{
    int n;
    cin>>n;
    while(n--)
    {
        int a,b,p;
        cin>>a>>b>>p;
        cout<<qmi(a,b,p)<<endl;
    }
    return 0;
}
//注意代码要放在两组三个点之间，才可以正确显示代码高亮哦~

作者：yankai
链接：https://www.acwing.com/activity/content/code/content/3748389/
来源：AcWing
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

（2）快速幂求逆元

推导：

$gif.latex?%5Cfrac%7Ba%7D%7Bb%7D%3Da*x%5C%20%28mod%5C%20m%29$

$gif.latex?%5Cfrac%7Ba%7D%7Bb%7D%3Da*b%5E%7B-1%7D%5C%20%28mod%5C%20m%29$

$gif.latex?b*b%5E%7B-1%7D%3D1%5C%20%28mod%5C%20m%29$

$gif.latex?bx%3D1%5C%20%28mod%5C%20m%29$

由费马小定理：

$gif.latex?b%5E%7Bp-1%7D%3D1%5C%20%28mod%5C%20p%29$

$gif.latex?b*b%5E%7Bp-2%7D%3D1%5C%20%28mod%5C%20p%29$

因此，x的一个可能解为：

b^p-2

求逆元，即用快速幂求出b的m-2次方模m的值。

代码：这里保证p是质数，因此可以用上述方法求解。如果不是质数，考虑扩展欧几里得算法。

//这里填你的代码^^
#include<iostream>
#include<algorithm>
using namespace std;

typedef long long LL;


LL qmi(int a,int b,int p)
{
    //考虑b的二进制形式
    LL res=1;
    while(b)
    {
        if(b&1) res=res*a%p;//如果b该位有1
        b>>=1;//去除一位
        a=a*(LL)a%p;//a平方
    }
    return res;
}


int main()
{
    int n;
    cin>>n;
    while(n--)
    {
        int a,p;
        cin>>a>>p;
        int res=qmi(a,p-2,p); 
        if(a%p) cout<<res<<endl;
        else puts("impossible");//不存在情况：a与p不互质，，模p必为0
        //ps：判断条件不能是res!=0，因为在p==2的时候会返回1
    }
    return 0;
}
//注意代码要放在两组三个点之间，才可以正确显示代码高亮哦~

作者：yankai
链接：https://www.acwing.com/activity/content/code/content/3748546/
来源：AcWing
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

二、扩展欧几里得算法

（1）推导和实现

扩展欧几里得算法：给定任意一对正整数a，b，求得一对整数x，y使得ax+by=gcd（a，b）。

推导：

回顾欧几里得算法，求最大公约数的边界时（即递归边界）：

$gif.latex?%28a%2Cb%29%3D%28b%2Ca%5C%20mod%5C%20b%29%3D.....%3D%28a%5E%7B%27%7D%2C0%29$

此时，对于（a‘，0），x，y的一组解为

$gif.latex?x%3D1%2Cy%3D0$

考虑完了边界条件，在考虑一步步返回x，y的方式。对于（b，a mod b）的y，x，与（a，b）的x，y

$gif.latex?y*b+%28a%5C%20mod%5C%20b%29*x%3Dd$

$gif.latex?b*y-%28a-%5B%5Cfrac%7Ba%7D%7Bb%7D%5D*b%29*x%3Dd$

$gif.latex?a*x+%28y-%5B%5Cfrac%7Ba%7D%7Bb%7D%5D*x%29*b$

从此可以看出，在一步步返回时，x不需要变动，而y需要按下式变化

$gif.latex?y%3Dy%5E%7B%27%7D-%5B%5Cfrac%7Ba%7D%7Bb%7D%5D*x$

拓展：实际上，x，y的解的结构为：

代码如下：

//这里填你的代码^^
#include<iostream>
using namespace std;

int egcd(int a,int b,int& x,int& y)
{
    if(!b)
    {
        x=1;
        y=0;
        return a;//如果b是0，a就是最大公约数，显然x=1，y=0是一组解
    }
    int d=egcd(b,a%b,y,x);//注意交换了y，x位置
    //此时yb+x(a%b)=d。。所以对于a，b。(y-a/b*x)b+ax=d;对于a，b来说需要变换一下y
    y -= a/b*x;
    return d;
}



int main()
{
    int n;
    cin>>n;
    while(n--)
    {
        int a,b,x,y;
        cin>>a>>b;
        egcd(a,b,x,y);
        cout<<x<<" "<<y<<endl;
    }
    return 0;
}
//注意代码要放在两组三个点之间，才可以正确显示代码高亮哦~

作者：yankai
链接：https://www.acwing.com/activity/content/code/content/3756046/
来源：AcWing
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

（2）线性同余方程

给定a，b，m求出一个x，使得 $gif.latex?a*x%3Db%5C%20mod%5C%20m$ 。这样的方程称为线性同余方程。

解法推导：

$gif.latex?a*x%3Db%5C%20mod%5C%20m$

$gif.latex?a*x+m*y%3Db$

由扩展欧几里得算法，我们可以求得k1，k2，使得：

$gif.latex?a*k_%7B1%7D+m*k_%7B2%7D%3Dgcd%28a%2Cm%29$

因此，由裴蜀定理知，b一定是gcd（a，m）的倍数，若不是则方程无解。

若有解，易知：

$gif.latex?x%3Dk_%7B1%7D*%5Cfrac%7Bb%7D%7Bgcd%28a%2Cm%29%7D$

代码：

//这里填你的代码^^
#include<iostream>
using namespace std;
typedef long long LL;

int egcd(int a,int b,int& x,int& y)
{
    if(!b)
    {
        x=1,y=0;
        return a;
    }
    int d=egcd(b,a%b,y,x);
    y-=a/b*x;
    return d;

}


int main()
{
    int n;
    cin>>n;
    while(n--)
    {
        int a,b,m;
        cin>>a>>b>>m;
        int x,y;
        int d=egcd(a,m,x,y);
        if(b%d) puts("impossible");
        else cout<<(LL)x*b/d%m<<endl;

    }
    return 0;
}
//注意代码要放在两组三个点之间，才可以正确显示代码高亮哦~

作者：yankai
链接：https://www.acwing.com/activity/content/code/content/3756228/
来源：AcWing
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

三、中国剩余定理

（1）定理阐述

中国剩余定理给定了以下的一元线性同余方程组：

$gif.latex?%28S%29%3A%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D%20%26%20x%5Cequiv%20a_%7B1%7D%5C%20%28mod%5C%20m_%7B1%7D%29%5C%5C%20%26%20x%5Cequiv%20a_%7B2%7D%5C%20%28mod%5C%20m_%7B2%7D%29%5C%5C%20%26%20x%5Cequiv%20a_%7B3%7D%5C%20%28mod%5C%20m_%7B3%7D%29%5C%5C%20%26.............................%20%5C%5C%20%26%20x%5Cequiv%20a_%7Bn%7D%5C%20%28mod%5C%20m_%7Bn%7D%29%20%5Cend%7Bmatrix%7D%5Cright.$

中国剩余定理说明：假设整数m1,m2, ... ,mn两两互质，则对任意的整数：a1,a2, ... ,an，方程组 $gif.latex?%28S%29$ 有解，且通解可以用如下方式构造得到：

设 $gif.latex?M%3D%5Cprod_%7Bi%3D1%7D%5E%7Bn%7Dm_%7Bi%7D$ ，并设 $gif.latex?M_%7Bi%7D%3D%5Cfrac%7BM%7D%7Bm_%7Bi%7D%7D%2C%5C%20%5Cforall%20i%5Cin%20%5Cbegin%7BBmatrix%7D1%2C2%2C3...n%20%5Cend%7BBmatrix%7D$ 。设 $gif.latex?t_%7Bi%7D%3DM%5E%7B-1%7D_%7Bi%7D$ ，为 $gif.latex?M_%7Bi%7D$ 模 $gif.latex?m_%7Bi%7D$ 的逆元。

则方程组通解为：

$gif.latex?x%3D%5Csum_%7Bi%3D1%7D%5E%7Bn%7Da_%7Bi%7Dt_%7Bi%7DM_%7Bi%7D+kM%2C%5C%20k%5Cin%20Z$

在模M的意义下，方程只有一个解：

$gif.latex?x%3D%5Csum_%7Bi%3D1%7D%5E%7Bn%7Da_%7Bi%7Dt_%7Bi%7DM_%7Bi%7D$

原理简述：m之间两两互质，因此模mi时，下标同的mj mod mi ==0，而mi*mi的逆元模mi为1。

（2）样题：

解题推导：

首先选出方程组中两个方程：

$gif.latex?%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D%20x%5C%20mod%20%5C%20a_%7B1%7D%5Cequiv%20m_%7B1%7D%5C%5C%20x%5C%20mod%20%5C%20a_%7B2%7D%5Cequiv%20m_%7B2%7D%20%5Cend%7Bmatrix%7D%5Cright.$

则：

$gif.latex?%5Cleft%5C%7B%5Cbegin%7Bmatrix%7Dx%3Dk_%7B1%7Da_%7B1%7D+m_%7B1%7D%5C%5C%20x%3Dk_%7B2%7Da_%7B2%7D+m_%7B2%7D%5Cend%7Bmatrix%7D%5Cright.$

$gif.latex?k_%7B1%7Da_%7B1%7D+m_%7B1%7D%3Dk_%7B2%7Da_%7B2%7D+m_%7B2%7D$

$gif.latex?k_%7B1%7Da_%7B1%7D+k_%7B2%7Da_%7B2%7D%3Dm_%7B2%7D-m_%7B1%7D$

可以利用扩展欧几里得算法求得该式中的一组k1，k2。根据线性同余方程组解的结构，k1和k2的通解为：

$gif.latex?%5Cleft%5C%7B%5Cbegin%7Bmatrix%7D%20K_%7B1%7D%3Dk_%7B1%7D+k*%5Cfrac%7Ba_%7B2%7D%7D%7Bd%7D%5C%5C%20K_%7B2%7D%3Dk_%7B2%7D+k*%5Cfrac%7Ba_%7B1%7D%7D%7Bd%7D%20%5Cend%7Bmatrix%7D%5Cright.%5C%20%5C%20%5C%20%5C%20%5C%20%5C%20%5C%20%5C%20k%5Cin%20Z$

$gif.latex?%5Ctherefore%20x%3Dk_%7B1%7Da_%7B1%7D+m_%7B1%7D%3D%28k_%7B1%7D+k%5Cfrac%7Ba_%7B2%7D%7D%7Bd%7D%29a_%7B1%7D+m_%7B1%7D$

$gif.latex?x%3Da_%7B1%7Dk_%7B1%7D+m_%7B1%7D+k%5Cfrac%7Ba_%7B1%7Da_%7B2%7D%7D%7Bd%7D%3Da_%7B1%7Dk_%7B1%7D+m_%7B1%7D+k%5Ba_%7B1%7Da_%7B2%7D%5D$

可看成：

$gif.latex?x%3D%28a_%7B1%7Dk_%7B1%7D+m_%7B1%7D%29+k%5Ba_%7B1%7Da_%7B2%7D%5D%3Dka+m_%7B0%7D$

这样就将两个方程合并为了一个方程，如此把所有方程合并起来，求得最后方程形式为：

$gif.latex?x%3Dka+m%20%5C%20%5C%20%5C%20k%5Cin%20Z$

此时令k合理取值，使得x为方程能够达到的最小正整数即位题目所求。

细节注意：

1.用扩展欧几里得算法求解时，注意k1要扩大一个相应的倍数，有可能无解。

2.过程中需要不断累积的k1可能爆int，即需要 $gif.latex?k_%7B1%7D+k%5Cfrac%7Ba_%7B2%7D%7D%7Bd%7D$ 方程需要控制k1为最小正整数解。

代码如下：

//这里填你的代码^^
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long LL;


LL egcd(LL a,LL b,LL &x,LL &y)
{
    if(!b)
    {
        x=1,y=0;
        return a;
    }
    int d=egcd(b,a%b,y,x);
    y-=a/b*x;
    return d;
}


int main()
{
    int n;
    cin>>n;
    LL a1,m1;
    cin>>a1>>m1;
    bool has_ans=true;

    for(int i=0;i<n-1;i++)
    {
        LL a2,m2;
        cin>>a2>>m2;
        LL k1,k2;
        LL d=egcd(a1,a2,k1,k2);
        if((m2-m1)%d) 
        {
            has_ans=false;
            break;
        }
        k1 *=(m2-m1)/d;//扩大倍数

        LL t=a2/d;
        k1=(k1%t+t)%t;//保留方程最小正整数解

        m1=a1*k1+m1;//合并方程
        a1=abs(a1/d*a2);

    }

    if(has_ans)
    {
        cout<<(m1%a1+a1)%a1;//保留方程最小正整数解
    }
    else cout<<"-1";
    return 0;
}
//注意代码要放在两组三个点之间，才可以正确显示代码高亮哦~

作者：yankai
链接：https://www.acwing.com/activity/content/code/content/3763166/
来源：AcWing
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。