zoukankan      html  css  js  c++  java
  • 拓展欧几里得算法

    先贴个模板: 
    1
    void gcd(LL a,LL b,LL &d,LL &x,LL &y) 2 { 3 if (!b) 4 { 5 x=1; 6 y=0; 7 d=a; 8 } 9 else 10 { 11 gcd(b,a%b,d,y,x); 12 y-=a/b*x; 13 } 14 }

     

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

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

    假设我们已经解出了 bx0 + (a%b)y0=gcd(a,b)   (1)

    令y=x0 -> ay0 + by = gcd(a,b)    (2)

    (2)-(1)得到 (a-a%b)y0 + b(y-x0)=0

    整理得到 y=x0-((a-a%b)/b)y0 = x0 - a/b*y0 = x0 - 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+x0) + b(y+y0)=d

    那么相减得到 ax0+by0=0  -> y0=-a/b * x0= -a'/b' * x0

    y0和x0都是整数,所以x0必定是b'的倍数. 即x0=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 m1=a1

    x mod m2=a2

    ...

    x mod mn=an

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

    设一个合法的解为t.

    那么t=m1*x+a1=m2*y+a2

    移项得m1*x+m2*y=a2-a1

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

    设m1*(x+x0)+a1=m2*(y+y0)+a2

    那么可以得到y0=m1/m2*x0

    设gcd(m1,m2)=d   m1'=m1/d   m2'=m2/d

    y0=m1'/m2' * x0

    所以x0必须是m2'的倍数。

    设之前解出来的x是x'

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

    那么所有符合要求的t=m* x'+a1 + k * m2' * m1 = t' + k * lcm(m1,m2).

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

    x mod lcm(m1,m2) = t'

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

    贴个模板:

     1 void equation()
     2 {
     3     cin >> m1 >> a1;
     4     for (int i=1;i<n;i++)
     5     {
     6         cin >> m2 >> a2; //t=m1*x+a1=m2*y+a2 -> m1*x+m2*y=a2-a1
     7         a=m1;
     8         b=m2;
     9         c=a2-a1;
    10         gcd(a,b,x,y,d);
    11         
    12         if (c%d!=0) //判断无解 
    13         {
    14            cout<<-1<<endl;
    15            return 0;
    16         }
    17         
    18         x*=c/d;  // 映射x 
    19         x%=m2/d;  // a1最终要mod lcm(m1,m2) t=m1*x+a1 x不能超过m2/d,否则容易爆long long 
    20         a1=m1*x+a1;
    21         m1=m1/d*m2;
    22         if (a1<0)  //令a1>=0 
    23            a1+=-a1/m1*m1+m1;
    24         if (a1>0)
    25            a1%=m1;
    26     }
    27 }
  • 相关阅读:
    iBatis的基本使用
    修改浏览器语言环境
    JavaScript实现级联下拉框
    Spring2.5与JDK8的集成问题
    JDK API文档下载
    C程序:年转化天
    线程基础
    拷贝控制之拷贝、赋值、销毁
    jquery中attr、prop、data
    input 类型为number型时,maxlength不生效?
  • 原文地址:https://www.cnblogs.com/vb4896/p/4009181.html
Copyright © 2011-2022 走看看