数论算法--快速幂、欧几里得算法、扩展欧几里得算法、求逆元、线性同余方程

快速幂

markdown中符号不太好打，举个例子：$3^{11}=3^{1011}=3^8+3^2+3^1$

代码如下：

typedef long long ll;
int qmi(int a,int k,int p){ // quickMi的缩写
    int res=1;
    while(k){
        if(k&1)res=(ll)a*res%p; // 乘法有可能溢出
        k>>=1;
        a=(ll)a*a%p;			// 有可能溢出
    }
    return res;
}

欧几里得算法（辗转相除法）

欧几里得算法的核心是$gcd(a,b)=gcd(b,a\%b)$

迭代代码：

int gcd(int a,int b){
    int t=a%b;
    while(t){
        a=b;
        b=t;
        t=a%b;
    }
    return b;
}

递归代码：

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

扩展欧几里得算法

裴蜀定理：如果$a,b$为正整数，则存在整数$x,y$使得$ax+by=gcd(a,b)$

一个例子：$gcd(252,198)=18，18=4\ast 252-5\ast 198$

迭代版本如下：
实例可看《离散数学及其应用（第八版）》例18

int exgcd(int a,int b,int &x,int &y){
    vector<int>s,f; 
    s.push_back(1);  // 插入的数是固定的
    s.push_back(0);
    f.push_back(0);
    f.push_back(1);
    vector<int>q;
    int t=a%b;
    q.push_back(a/b);
    while(t){
        a=b;
        b=t;
        t=a%b;
        q.push_back(a/b); 
    }  // 求最大公约数
    for(int i=2;i<=q.size();i++){
        int s_temp=s[i-2]-s[i-1]*q[i-2];
        s.push_back(s_temp);
        int t_temp=f[i-2]-f[i-1]*q[i-2];
        f.push_back(t_temp);
    }
    x=s.back(),y=f.back();
}

递归版本：

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

逆元

除法比较复杂，最好转化为除法，若$\frac{a}{b}\equiv ax(modm)$，则$x$为$b$的逆元。
上式可变形为变形为 $bx\equiv 1(modm)$，

扩展欧几里得求逆元

继续变形 $bx-my\equiv 1（modm）$，转化为 $bx+m{y}^{{}^{\prime }}\equiv 1(modp)$，这个可以用扩展欧几里得算法 $ax+by\equiv gcd(a,b)$，来求解，只不过这里 $gcd(a,b)=1$。

int exgcd(int a,int b,int &x,int &y){ // 扩展欧几里得算法
    if(!b){
        x=1,y=0;
        return a;
    }
    int d=exgcd(b,a%b,y,x);
    y-=(ll)a/b*x;
    return d;
}

void mod_reverse(int a,int m){  	// 求逆元
    int x,y;
    int d=exgcd(a,m,x,y);
    if(d!=1)cout<<"不能构成逆元"<<endl;  // 判断a与m的最小公倍数是否为1,若为1，则有逆元，若不为1，则没有逆元。
    else cout<<"逆元为"<<(x%m+m)%m<<endl;
}

快速幂求逆元

快速幂求逆元需要由一定的限制条件。
$bx\equiv 1(modm)$
由费马小定理可知：当$b$为质数时，$b^{n-1}\equiv 1(modm)$
变形为$b\ast b^{n-2}\equiv 1(modm)$，可得$x=b^{n-2}$，即b的逆元为$b^{n-2}$（$b,m$互质，即$a\%p!=0$）

代码如下：

int qmi(int a,int k,int p){ // 求快速幂
    int res=1;
    while(k){
        if(k&1)res=(LL)res*a%p;
        k>>=1;
        a=(LL)a*a%p;
    }
    return res;
}

void mod_reverse(int a,int p){  // 求逆元
    if(a%p){	// a,p互质
        cout<<"逆元为:"<<qmi(a,p-2,p)<<endl;
    }
    else cout<<"不存在逆元"<<endl;
}

线性同余方程

形如$axb\equiv (modm)$，其中$m$为正整数，$a,b$为整数，$x$为变量。这样的方程，称为线性同于方程。
变形$ax\equiv my+b$
–>$ax-my\equiv b$
–>$ax+m{y}^{{}^{\prime }}=b$
又因为$ak+ml=gcd(a,m)$，因此$b$一定是$gcd(a,b)$的倍数，因为$a$是$gcd(a,m)$，$m$也是$gcd(a,b)$的倍数，所以他们的倍数相加$ax+my$也是$gcd(a,b)$的倍数。可通过扩展欧几里得算法求。最后判断$b\%a==0$，然后乘$b/gcd(a,m)$即可。

代码如下：

int exgcd(int a,int b,int &x,int &y){  // 扩展欧几里得算法
    if(!b){
        x=1,y=0;
        return a;
    }
    int d=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return d;
}

void linear(int a,int b,int m){  //解线性同余方程
    int x,y;
    int d=exgcd(a,m,x,y);
    if(b%d)cout<<"不存在上述方程"<<endl;
    else cout<<(ll)x*(b/d)%m<<endl; //乘 b/d。
}