初等数论 ————拓展欧几里得算法

拓展欧几里得（exgcd）简直在数论中无比重要的概念‘

上文说到gcd可以在时间复杂度为log（2 n）下求出两个数的最大公因数

exgcd我知道就两个 求二元一次方程方程（也可以叫不定方程）的通解  和逆元

贝祖定理

贝祖定理是代数几何中一个定理，其内容是若设a,b是整数，则存在整数x,y，使得ax+by=gcd（a,b），(a,b)代表最大公因数，则设a,b是不全为零的整数，则存在整数x,y，使得ax+by=(a，b)。

————————————————————————————————————————————————

对于不完全为 0 的非负整数 a，b，gcd（a，b）表示 a，b 的最大公约数，必然

存在整数对 x，y ，使得 gcd（a，b）=ax+by

求解 x，y的方法的理解

设 a>b。

1，显然当 b=0，gcd（a，b）=a。此时 x=1，y=0；

2，a>b>0 时

设 ax1+ by1= gcd(a,b);

bx2+ (a mod b)y2= gcd(b,a mod b);

根据朴素的欧几里德原理有 gcd(a,b) = gcd(b,a mod b);

则:ax1+ by1= bx2+ (a mod b)y2;

即:ax1+ by1= bx2+ (a - [a / b] * b)y2=ay2+ bx2- [a / b] * by2;

说明： a-[a/b]*b即为mod运算。[a/b]代表取小于a/b的最大整数。

也就是ax1+ by1 == ay2+ b(x2- [a / b] *y2);

根据恒等定理得：x1=y2; y1=x2- [a / b] *y2;

这样我们就得到了求解 x1,y1 的方法：x1，y1 的值基于 x2，y2.

上面的思想是以递归定义的，因为 gcd 不断的递归求解一定会有个时候 b=0，所以递归可以结束

举个例子

已知不定方程为

  

，利用辗转相除法求出一组整数解 

解：求  的算式为：

所以

即

所以

所以

  

是不定方程

  

的一组解

#include<iostream>
using namespace std;
int exgcd(int a,int b,int &x,int &y)
{
    if(b==0)
    {
        x=1;y=0;
        return a;
        cout<<"*******"<<endl;
    }
    int r=exgcd(b,a%b,x,y);
    int t=x;x=y;y=t-a/b*y;
    cout<<x<<" "<<y<<endl;
    //cout<<"#############"<<endl;
    return r;
}
int main()
{
    int a,b,x,y,r;
    cin>>a>>b>>x>>y;
    r=exgcd(a,b,x,y);
    cout<<r<<endl;
    return 0;
}

上面已经列出找一个整数解的方法，在找到p * a+q * b = Gcd(a, b)的一组解p0,q0后，p * a+q * b = Gcd(a, b)的其他整数解满足：

p = p0 + b/Gcd(a, b) * t

q = q0 - a/Gcd(a, b) * t(其中t为任意整数)     //这里我不明白

至于pa+qb=c的整数解，只需将p * a+q * b = Gcd(a, b)的每个解乘上 c/Gcd(a, b) 即可，但是所得解并不是该方程的所有解，找其所有解的方法如下：

在找到p * a+q * b = Gcd(a, b)的一组解p0,q0后，可以

得到p * a+q * b = c的一组解p1 = p0*(c/Gcd(a,b)),q1 = q0*(c/Gcd(a,b))，p * a+q * b = c的其他整数解满足：

p = p1 + b/Gcd(a, b) * t

q = q1 - a/Gcd(a, b) * t(其中t为任意整数)

p 、q就是p * a+q * b = c的所有整数解。

 逆元（Multiplicative inverse modulo）

如果ax≡1 (mod p),且gcd(a,p)=1（a与p互质），则称a关于模p的乘法逆元为x

因为取模不满足除法 所以定义在mod p的环境中 某个数除以a等于乘以x

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

设c是b的逆元，则有b*c≡1(mod m)；

则(a/b)%m = (a/b)*1%m = (a/b)*b*c%m = a*c(mod m);

即a/b的模等于a*b的逆元的模；

#include <iostream>
#include <cstdio>
using namespace std;
int x,y,q;
void extend_Eulid(int a,int b)
{
if(b == 0){
x = 1;y = 0;q = a;
}else{
extend_Eulid(b,a%b);
int temp = x;
x = y;
y = temp - a/b*y;
}
}
int main()
{
int a,b;
cin>>a>>b;
extend_Eulid(a,b);
x = (x % p + p) % p;//如果x为负数则需要这一步 
printf("%d=(%d)*%d+(%d)*%d
",q,x,a,y,b);
return 0;
}

 给定 a 和b。
a 要有逆元 ， 那么gcd（ a ， b ) = 1
假设a的逆元 为x ， 那么就有 a * x mod b = 1
也就是 a * x + b * y = 1
上面的程序中输入a和b就可以求出对应的x和y。
其中 x 就是 a的逆元

 贴个例题http://poj.org/problem?id=1061

#include <iostream>
#define dnt long long
using namespace std;
dnt x, y;
dnt a, b, m, n, L;
dnt Exgcd( dnt a, dnt b, dnt &x, dnt &y ) {
    if ( b == 0 ) {
        x = 1;
        y = 0;
        return a;
    }
    dnt d =  Exgcd(b, a%b, x, y), temp = x;
    x = y;
    y = temp-a/b*y;
    return d;
}

dnt solv( dnt a, dnt b, dnt c ) {
    dnt d = Exgcd(a, b, x, y);
    if ( c % d ) return -1;
    x = x * c / d;
    y = y * c / d;
    x = (x % b + b) % b;
    return x;
}
int main() {
    cin >> a >> b >> m >> n >> L;
    if ( solv(n-m, L, a-b) != -1 )
    cout << solv(n-m, L, a-b) << endl;
    else cout << "Impossible" << endl;
    return 0; 
}