zoukankan      html  css  js  c++  java
  • 模线性方程(递归版+迭代版)& 扩展欧几里德

      线性方程:设ab是两个整数,g = gcd(a,b)ab的最大公约数。求满足方程 a*x + b*y = g 的整数解xy

      递归版:扩张欧几里德

      在用欧几里德算法算ab的最大公约数时,我们依次得到:

      a = q(1) * b + r(1)

      b = q(2) * r(1) + r(2)

      r1 = q(3) * r(2) + r(3)

        ......

      r(n-2) = q(n) * r(n-1) + r(n)

      当r(n)0r(n-1)便是ab的最大公约数g。我们在Euclid递归返回的过程中,构造解。

      设第k-1层有:

        r(k-3) = q(k-1) * r(k-2) + r(k-1)

      第k层有:

        r(k-2) = q(k) * r(k-1) + r(k)

      递归时,设第k层的解已经构造好了,即我们已经知道了:x(k)y(k),满足:x(k) * r(k-2) + y(k) * r(k-1) = g.

      则第k-1层为:x(k-1) * r(k-3) + y(k-1) * r(k-2) = g,即x(k-1)*( q(k-1)*r(k-2) + r(k-1) ) + y(k-1)*r(k-2) = g。整理得:x(k-1)*r(k-1) + (x(k-1)*q(k-1) + y(k-1))*r(k-2) = g

      得到系数的转移方程:

        x(k-1) = y(k)

        y(k-1) = x(k) - y(k)*q(k-1)

      扩展欧几里德代码:

     1 //zzy2012.7.17
     2 #include<cstdio>
     3 #include<iostream>
     4 
     5 using namespace std;
     6 
     7 int Extended_Euclid(int a,int b,int &x,int &y){
     8     int q,temp,g;
     9     if(b==0){
    10         x = 1;
    11         y = 0;
    12         return a;
    13     }
    14     g = Extended_Euclid(b,a%b,x,y);
    15     q = a/b;
    16     temp = x;
    17     x = y;
    18     y = temp - q*y;
    19     return g;
    20 
    21 }
    22 
    23 int main()
    24 {
    25     int a,b,x,y,g;
    26     cin>>a>>b;
    27     g = Extended_Euclid(a,b,x,y);
    28     cout<<a<<'*'<<x<<'+'<<b<<'*'<<y<<'='<<g<<endl;
    29     return 0;
    30 }

      迭代版

      考察第k-1层:

        r(k-3) = q(k-1) * r(k-2) + r(k-1)

      得r(k-1) = r(k-3) - q(k-1) * r(k-2),得系数递推关系:

        x(k) = x(k-2) - q(k) * x(k-1)

        y(k) = y(k-2) - q(k) * y(k-1)

      迭代版代码:

     1 //zzy2012.2.21
     2 //求解形如a*x+b*y=gcd(a,b)的线性方程
     3 #include<cstdio>
     4 #include<iostream>
     5 
     6 using namespace std;
     7 
     8 typedef struct
     9 {
    10     int x,y;
    11 }node;
    12 
    13 node r[3];
    14 
    15 int solve(int a,int b)
    16 {
    17     int f,q,t;//f指向r元素的下标,q每次迭代算出的商,t保存余数,r[]数组保存每次计算出的系数
    18     f=2;
    19     r[0].x=1;
    20     r[0].y=0;
    21     r[1].x=0;
    22     r[1].y=1;
    23     while(b!=0)
    24     {
    25         q=a/b;
    26         t=a%b;
    27         a=b;
    28         b=t;
    29         r[f%3].x=r[(f-2)%3].x-q*r[(f-1)%3].x;
    30         r[f%3].y=r[(f-2)%3].y-q*r[(f-1)%3].y;
    31         f++;
    32     }
    33     return (f-2)%3;
    34 }
    35 
    36 int main()
    37 {
    38     int a,b,f;
    39     cin>>a>>b;
    40     if(a!=0 || b!=0)
    41     {
    42         f=solve(a,b);
    43         printf("%d %d\n",r[f].x,r[f].y);
    44     }
    45     return 0;
    46 }

       同余式a*x = c (mod m)可以转换为求线性方程。

      代码:

     1 /*-----ЭЌгрЪН-----2011.7.10-----Library-----zzy-----*/
     2 #include<cstdio>
     3 #include<cmath>
     4 
     5 int arg[2];
     6 
     7 void ini(int a,int b)
     8 {
     9     arg[0]=0;
    10     arg[1]=1;
    11     arg[2]=1;
    12     arg[3]=-a/b;
    13 }
    14 
    15 int lineeq(int a,int b)
    16 {
    17     int x,q;
    18     if(b==0)
    19         return a;
    20     q=-a/b;
    21     x=arg[0]-q*arg[1];
    22     arg[0]=arg[1];
    23     arg[1]=x;
    24     return lineeq(b,a%b);
    25 }
    26 
    27 int main()
    28 {
    29     int a,c,m,x,g,i;
    30     scanf("%d %d %d",&a,&c,&m);
    31     if(m==0)
    32     {
    33         printf("illegal m.\n");
    34         return 0;
    35     }
    36     if(a==0)
    37     {
    38         if(c%m==0)
    39             printf("x can be any integer.\n");
    40         else
    41             printf("no solution.\n");
    42         return 0;
    43     }
    44     ini(a,-m);
    45     g=lineeq(-m,a%(-m));
    46     g=abs(g);
    47     if(c%g!=0)
    48     {
    49         printf("no solution.\n");
    50         return 0;
    51     }
    52     x=arg[0]*c/g;
    53     for(i=0; i<g; i++)
    54     {
    55         printf("%d + k * %d\n",x+i*m/g,m);
    56     }
    57     return 0;
    58 }
  • 相关阅读:
    第三次上机练习
    第三次作业
    第二次上级练习
    第二次作业
    第一次上机练习
    第一次作业
    4.20
    4.16
    4.10
    4.9
  • 原文地址:https://www.cnblogs.com/Lattexiaoyu/p/2595782.html
Copyright © 2011-2022 走看看