拓展欧几里得算法

先贴个模板： 
1 void gcd(LL a,LL b,LL &d,LL &x,LL &y)
{
    if (!b)
    {
        x=1;
        y=0;
        d=a;
    }
    else
    {
        gcd(b,a%b,d,y,x); 
        y-=a/b*x;
    }
}

 

下面给出我对这个算法的理解:

现在要解ax+by=gcd(a,b)

假设我们已经解出了 bx₀ + (a%b)y₀=gcd(a,b)   (1)

令y=x₀ -> ay₀ + by = gcd(a,b)    (2)

(2)-(1)得到 (a-a%b)y₀ + b(y-x₀)=0

整理得到 y=x₀-((a-a%b)/b)y₀ = x₀ - a/b*y₀ = x₀ - a/b * x

根据上面的递归,返回的时候x=y0,y=x0 所以y需要再减去 a/b*x

函数结束时就可以得到一组解了。那么如何得到其他解呢？

首先令 d=gcd(a,b)   a'=a/d  b'=b/d

ax+by=d

设a(x+x₀) + b(y+y₀)=d

那么相减得到 ax₀+by₀=0  -> y₀=-a/b * x₀= -a'/b' * x₀

y₀和x₀都是整数,所以x₀必定是b'的倍数. 即x₀=k*b'

因此对于方程ax+by=d , 让x=x+b' y=y-a' 方程仍然成立。

所以方程的解为

x=x+k*b'

y=y- k*a' 

 

 

为了便于理解，如果给出的a<0,一开始就把a=-a，c=-c，这样a，b，c都是正的了。

如果 （c%d！=0） 那么无解。

否则 令k=c/d， a'=a/d ,   b'=b/d;

对于函数求出的x，先映射一下 ，x*=k.

然后 x每次 加 b', y每次减 a' 是不会影响等式成立的。

这样 就可以 求出x的最小整数解。

     if (x<0)
        (x+=abs(x)/b' * b' + b') %= b';
     if (x>0)
        x-=x/b' * b';

 

扩展欧几里得还可以拿来求解同余方程组:

x mod m₁=a₁

x mod m₂=a₂

...

x mod m_n=a_n

 

可以2个两个来合并,先看前2个：

设一个合法的解为t.

那么t=m₁*x+a₁=m₂*y+a₂

移项得m₁*x+m₂*y=a₂-a₁

用扩展欧几里得可以得到一组合法的x,y

设m₁*(x+x₀)+a₁=m₂*(y+y₀)+a₂

那么可以得到y₀=m₁/m₂*x₀

设gcd(m₁,m₂)=d   m₁'=m₁/d   m₂'=m₂/d

y₀=m₁'/m₂' * x₀

所以x₀必须是m₂'的倍数。

设之前解出来的x是x'

那么所有符合要求的x=x' + k*m₂' 把x'带入得到一个合法的t'

那么所有符合要求的t=m₁ * x'+a₁ + k * m₂' * m₁ = t' + k * lcm(m₁,m₂).

因此2个同余方程就可以合并成一个了：

x mod lcm(m₁,m₂) = t'

两两合并就可以解出来了。

 

贴个模板：

 

void equation()
{
    cin >> m1 >> a1;
    for (int i=1;i<n;i++)
    {
        cin >> m2 >> a2; //t=m1*x+a1=m2*y+a2 -> m1*x+m2*y=a2-a1
        a=m1;
        b=m2;
        c=a2-a1;
        gcd(a,b,x,y,d);
        
        if (c%d!=0) //判断无解 
        {
           cout<<-1<<endl;
           return 0;
        }
        
        x*=c/d;  // 映射x 
        x%=m2/d;  // a1最终要mod lcm(m1,m2) t=m1*x+a1 x不能超过m2/d,否则容易爆long long 
        a1=m1*x+a1;
        m1=m1/d*m2;
        if (a1<0)  //令a1>=0 
           a1+=-a1/m1*m1+m1;
        if (a1>0)
           a1%=m1;
    }
}