zoukankan      html  css  js  c++  java
  • [数学]求解一类方程的方法 博文笔记&&扩展欧几里得、中国剩余定理学习笔记。

    原文传送:https://www.luogu.com.cn/blog/1239004072Angel/post-shuo-xue-qiu-xie-yi-lei-fang-cheng-di-fang-fa

      不推荐只看原文进行学习。

    裴蜀定理 (ax+by的最小正值还能整除a,b, 该最小正值为gcd(a,b))。

     方程ax+by=c (a,b,c为整数)有整数解的充要条件为gcd(a,b) | c 原博客只证明了充分性,这里证明必要性:

      若方程有整数解x1,y1,则c=ax1+by1。因为ax1+by1能被 gcd(a,b)整除,则c必然能被gcd(a,b)整除,否则方程无解。(一个方程变形罢了)

    扩展欧几里得算法:

      总的来说,扩展欧几里得用来求不定方程ax+by=c的整数解或线性同余方程ax≡c mod b。利用欧几里得算法框架递归到b=0后返回得到答案。

        推到a′x+b′y=b′x'+(a′modb′)y' ① 时(注意等号两边的a'是同一个数(b'同理),但x和x'不是同一个数(y和y'同理),原博客没有体现出来),可令a''=b',b''=a′modb′,那么方程能像这样再一直等于下去,直到 某个b'...'  (这里省略号是若干'的意思)为0 (a,b的变化就是欧几里得算法的过程),那么此时方程为a'...' x = gcd(a'...',0)=a'...',得到x'...'=1,y'...'=任意整数 的解。为了方便,可以取y'...'=0。加上 下文已知, 后面的方程的解可推出前面方程的解,那么一开始的方程就解出来了。

        一个更清楚的口胡解释:方程a′x+b′y=gcd(a',b')和b′x'+(a′modb′)y'=gcd(b′,a′modb′) 都有解(裴蜀定理)(也许用表示联立的大括号表示更直白),同时 gcd(a',b')=gcd(b′,a′modb′) 。解出方程b′x'+(a′modb′)y'=gcd(b′,a′modb′) 后,方程a′x+b′y=gcd(a',b')=gcd(b′,a′modb′)=b′x'+(a′modb′)y'  (此时x' 、y‘ 确定,b′x'+(a′modb′)y'为一个定值),这时方程最右边可变形为:a′y’+b′(x‘−b′/a′  ​ * y’),即​a′x+b′y=a′y’+b′(x‘−b′/a′  ​ * y’),这里面a',b',x',y'都已知,那么取x=y',y=x‘−  b′/a′  ​* y’就是原方程的一个解。这时可以知道,要求当前不定方程ax+by=gcd(a,b)的解,可看它对应的下一个方程bx'+(a%b)y'=gcd(b,a%b)的解,而这”下一个方程“的解,又可以看它的下一个方程即“下下个方程‘’…直到最后看的方程的第二个系数为零,得解,再按顺序返回到第一个方程,就得到了答案。(递归的依据就是gcd(a,b)=gcd(b,a%b),故叫扩展欧几里得,还能顺便求出最大公因数)

      

      对于方程ax+by=c  (gcd(a,b)|c) 设k=c/gcd(a,b)  求出ax+by=gcd(a,b)的解x,y后将x,y都*k 即为ax+by=c 的一组特解。

      代码:

     1 int exgcd(int a,int b,int &x,int &y)
     2 //调用函数完成后传址&x,&y的变量的值就为关于当前函数a,b的值的ax+by=gcd(a,b)的解(x,y) 
     3 {
     4     if(!b)
     5     {
     6         x=1;y=0;
     7         return a;
     8     }
     9     int g=exgcd(b,a%b,y,x);//反着传址,省去后文交换x,y的值 
    10     y-=a/b*x;
    11     return g;
    12 }
    13 
    14 int main()//求解不定方程ax+by=c
    15 {
    16     int a=read(),b=read(),c=read(),x,y;
    17     int g=exgcd(a,b,x,y),e=c%g;//按址传递形参。
    18     //不写按址传递形参的话,可以写结构体或全局变量  
    19     if(e)//无解 
    20     {
    21         cout<<"NONE"; 
    22         return 0; 
    23     }
    24     x*=c/g;y*=c/g;//不要忘了c不一定是gcd(a,b) 
    25     cout<<x<<endl<<y;
    26     return 0;
    27 }

        上面讨论的情况只是a,b,c都大于0的情况。当a,b,c有小于0的情况时,特判一下:若是a或b小于0,则完全可以把负号去掉,最后再给x或y变下号就好;若是c小于0,那么c/g小于0,并不需要做什么多余操作;若a,b小于0,可等号两边都提出一个-1来,转化成c小于0的情况;若a,b,c都小于0,则可等号两边都提出一个-1来,转化成a,b,c大于0都情况。第一种情况相当于换元,故x或y最后要额外处理,后两种情况是方程的变形,不会改变解。

        也可以不用这么麻烦。若会小于0,算法向下递归的正确性不变,只是递归到b=0时,a可能为负值,此时2种处理:

          1、算法返回a的绝对值,满足exgcd返回值仍为gcd(a,b),令x=  a<0 ?-1 :1   (因为gcd一定是非负的)

          2、令x=1,算法返回a,exgcd返回值(设为g)的绝对值才是gcd,且此时x,y是满足ax+by=g的一组解(ax+by=-gcd一样有解,那么也可凭此递归下去,因为能递归下一层的依据就是构造了两个等号右边一样(不管是gcd还是-gcd),左边系数有关系、可变形为同样形式的方程,这样求出一个方程的解,另一个就可根据形式写出解)。

      求出不定方程ax+by=c (a,b>0)的一组特解x0,y0后,他的所有解都可以表示为:x=x0+k*(b/gcd(a,b)),y=y0-k*(a/gcd(a,b)),k为整数。对于a,b有小于0的情况,根据上文来转化为a,b大于0的情况,最后决定是否对x,y变号。

        简单证明:不定方程ax+by=c (a,b>0)的解相对于x0,y0都可以表示成x=x0+u,y=y0+v的形式,且u和v异号(因为a,b>0)。设d=a,b的最小公倍数,则d=a*b / gcd(a,b)  。
          考虑x0缓慢增加(减少同理)以找到新的解,发现只有ax增加的量w为d的倍数时,by才能没有在b、y都为整数的情况下抵消ax的增加量(w一定为a的倍数,否则x不为整数,w若不为b的倍数,则y不为整数。故若要求解为整数(不定方程的要求),w必为a,b的公倍数,即为d的倍数,w=kd,同时也说明相对于x0,y0的ax的变化量不为a,b的公倍数的解x,y是不存在的)ax变化了kd(增加还是减少要看k的符号)=k*  a*b / gcd(a,b),则x变化了k* b/gcd(a,b),by要向反方向变化kd ,得y变化了 -k * a/gcd(a,b) 。得证x=x0+k*(b/gcd(a,b)),y=y0-k*(a/gcd(a,b))一定还是原方程的解

      

      时间复杂度:与欧几里得算法同级,上限都为logn,随机数据下会更快。

        简单证明:考虑a,b,a%b,b%(a%b)。a%b一定<b,猜想复杂度logn的话,考虑b%(a%b)与b/2的大小关系,而这又与a%b有关(因为b%(a%b) < a%b)。因为关键节点就是b/2,想到分两种情况:

          a%b>  b/2,则b%(a%b)=b- a%b < b/2;

          a%b<  b/2,则b%(a%b)<a%b< b/2;

        发现gcd或exgcd函数每递归两层,函数的第二参数的规模就减少至少一半,故复杂度上限级别为log b。

    例题:洛谷P2421 [NOI2002] 荒岛野人 https://www.luogu.com.cn/problem/P2421

      注意当ax+by=c的(a,b)=1时,通解x=x0+k*b的形式才成立,否则(a,b)不一定为1时通解应为标准形式x=x0+k*b/gcd(a,b),不要忘了一般情况下的通解公式的后部是有除以gcd的(要找全解,ax就要不断变化a,b最小公倍数的量(a,b才是常数),x就变化b/gcd(a,b)的量)。y同理。

      还是容易在做完exgcd后忘记给解乘c/gcd(a,b),别忘了做完exgcd求出的是方程ax+by=gcd(a,b)的解,是不是ax+by=c的解还要再详看(有可能还无解呢)。

      所谓的“想当然”最好自证一下再用。相似的概念并行思考时容易混淆,注意区分。

      注意:若想防止溢出,exgcd中不能简单地对x,y模b/d,a/d。通解公式中x,y一个变化了,另一个也会变化。  

    上文‘%’号未经特殊说明的话都指取余运算(毕竟我的语言环境是c++)。

    中国剩余定理CRT(本文该模块中bi是同余号右边的数,ai是模数,ti是mi在模ai下的逆元):

      总的来说,用于求解模数互质系数为1的线性同余方程组。利用各模数与相应逆元构造了一个某模数下对应项结果为bi,其他项全是0的可行解。

      原文的“数学归纳”原来也是假设最后的解有多个再相减证明公倍数整除差值。

      证明及详解(更清楚明确,还有个明确的性质):孙子定理  https://baike.baidu.com/item/%E5%AD%99%E5%AD%90%E5%AE%9A%E7%90%86/2841597?fromtitle=%E4%B8%AD%E5%9B%BD%E5%89%A9%E4%BD%99%E5%AE%9A%E7%90%86&fromid=11200132&fr=aladdin

      

       对于:   

        能在biti处模ai的原因:不管怎样,mi必须要是除了自己(ai)的其他所有模数的乘积,以保证在被其他模数取模时结果为0,只剩下那个模数对应的项,才能保证正确性。而biti在模ai时才会看到,在最后答案中模ai的话,不会影响正确性,同时减少数值范围(但不能对mi取模)

      

      代码很简朴,先求出模数积(O(n)),依次求出所有mi  (O(n)),ti (O(nlogn)),最后求出答案。总时间复杂度:O(nlogn)

      若方程组其中的同余方程的x前有系数,可用扩展欧几里得解出该方程,转化为系数为1的形式,最后总时间复杂度不变。

    扩展中国剩余定理(模数不一定两两互质时求解同余方程组):

      用于求解模数不一定互质的线性同余方程组。通过联立逐步合并方程组至一个方程。

      为方便理解中国剩余定理,有一个引理:对于两个互质的正整数模数 n, m,给出两个整数 a,b满足 0≤a<n,0≤b<m,则在[0,nm)中存在且仅存在一个 同余方程组{ x modn =a,x mod m=b }的解 。

        证明很容易理解,n,m互质时同余方程组一定有多个解(由中国剩余定理百度百科可知)。对于同余方程组的两个解x1,x2。作差得y ,y可同时被n 、m整除,即y为n,m的公倍数,那么在x1,x2他们的差距一定是n,m最小公倍数的倍数,x1,x2不相等的情况下最少也差lcm(n,m)(n,m互质的情况下lcm(n,m)即为nm),此时引理结论可得证。

        若模数不一定互质(对应扩展中国剩余定理),在方程组有解时(有解的条件下文再讨论),也是会有多个解,对于这时的解 x1,x2,差y可同时被n 、m整除,仍为n,m的公倍数,x1,x2不相等的情况下最小差距还是n,m的最小公倍数。最后可知模数不一定互质时,若有解,则[0,lcm(n,m))中只会存在一个解 。

        由于 ai​ 之间可能不互质,故我们可能找不到某个 mi​ 在模 ai​ 意义下的逆元,就无法用CRT的公式求解。

      算法核心:变形等式,合并方程:

        

         两个方程组是等价的,对第二个方程组找条件性质/变形 求解,即为原方程组的解。

        又等价于:

         

         方程通解可由上文引理印证,也知lcm是通解的最小变化量。其实上图第一行变形的等号右边应该是b2-b1,代码中写的也是b2-b1。exgcd中的参数可以直接传a1和a2(不用传-a2),这样得到的k1,k2乘上(b2-b1)/exgcd 后,k1是正确值,k2是正确值的相反数,但要求出x,只需求出k1和k2中的一个正确值就好。

        故对于 n 个方程,我们顺次合并两个方程即可。在合并过程中发现一处无解,那么方程组就无解。

      时间复杂度:nlogn

      启发思路:

        对于有很多个方程的方程组,可以划分部分,在各个部分里求出满足当前部分的所有方程的解,最后再求同时满足各部分的解的解,即为总方程组的解(分治再合并);

        以后在整数意义下有两个未知数的方程都可解了。

      杂想:扩展CRT会取代中国剩余定理吗?在OI中,扩展CRT更泛用,运算时的数值范围大头在求lcm上,比普通CRT不需要实现判断模数两两互质、更不易溢出数据类型。但要明白,能逐个合并方程的前提是我们用计算机写代码,让计算机做这些重复的活,而中国剩余定理直接给出一个最后答案的表达式,适合人工计算,这也是日常中现阶段数学和信息的一个区别吧。

      例题:P4777 【模板】扩展中国剩余定理(EXCRT)https://www.luogu.com.cn/problem/P4777

      AC代码:

     1 #include<iostream>
     2 #include<cstdio>
     3 #include<cstring>
     4 #include<cmath>
     5 
     6 #define ll long long
     7 #define ull unsigned long long
     8 #define min(a,b) ((a)>(b)?(b):(a))
     9 #define max(a,b) ((a)>(b)>(a):(b))
    10 #define swap(a,b) ((a)^=(b),(b)^=(a),(a)^=(b))
    11 
    12 using namespace std;
    13 
    14 const int N=1e5+5;
    15 
    16 ll n,ai[N],bi[N],a1,b1,m,k1,k2; 
    17 
    18 inline ll read()//int N+
    19 {
    20     ll x=0;
    21     char ch=getchar();
    22     while(!isdigit(ch)) ch=getchar();
    23     while(isdigit(ch)) x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
    24     return x;
    25 }
    26 
    27 ll exgcd(ll a,ll b,ll &x,ll &y)
    28 {
    29     if(!b)
    30     {
    31         x=1;y=0;//此题a,b都非负 
    32         return a;
    33     }
    34     ll d=exgcd(b,a%b,y,x);
    35     y-=a/b*x;
    36     return d;
    37 }
    38 
    39 inline ll mul(ll a,ll b,ll mod)//快速乘,防溢出
    40 {
    41     ll ret=0;
    42     bool f=0;
    43     if(b<0)
    44     {
    45         f=1;
    46         b=-b;
    47     }
    48     while(b)
    49     {
    50         if(b&1)
    51             ret=(ret+a)%mod;
    52         b>>=1;
    53         a=(a+a)%mod;
    54     }
    55     return f?-ret:ret;
    56 }
    57 
    58 int main()
    59 {
    60     n=read();
    61     for(int i=1;i<=n;++i)
    62         ai[i]=read(),bi[i]=read();
    63     m=a1=ai[1],b1=bi[1];
    64     ll c,d;
    65     for(int i=2;i<=n;++i)//EXCRT算法核心过程 
    66     {
    67         c=bi[i]-b1;
    68         d=exgcd(a1,ai[i],k1,k2);
    69         m*=ai[i]/d;
    70         b1=(mul(mul(k1,c/d,m),a1,m)+b1)%m;//(k1*(c/d)*a1+b1)%m
    71         a1=m;
    72     }
    73     if(b1<0)
    74         b1+=m;//答案要>=0 
    75     printf("%lld",b1);
    76     return 0;
    77 }

        思路很简单,对溢出的处理废了点思考时间。说实话该代码没有对exgcd中做防溢出处理,luogu的数据还是弱了。除了高精度(极限情况下数值范围是个2*18*log1018 =2176位数,)有正确性更优的算法:https://www.luogu.com.cn/blog/user31435/solution-p4777

        思考过程想到的性质:

          lcm(a,b) >=lcm(b,a%b) 当a<b时等号成立

          若想防止溢出,exgcd中不能简单地对x,y模b/d,a/d。通解公式中x,y一个变化了,另一个也会变化。

    杂项:

      两算法都不要忘判断无解的情况。扩展欧几里得初步求出的解的是ax+by=gcd(a,b)的解,不代表ax+by=c一定有解。

      做数据范围很大的乘法时可能要注意模不同数时不同的应用环境,例如在CRT解的公式的最终结果是对M取模,除了前文提到的,计算时不能随便对某个模数处处取模。很多式子后面没有注明“模的意义下”,不要一看到大数计算就想当然要取模,会碍到/误导 思考分析。

      为什么在取模意义下的乘法式子中无论以任意顺序取模,结果都是正确且唯一的:设模数为m,一个数/式子的结果都能分解成km+r的形式,由分配律拆开后,km的那部分无论乘什么,结果在模意义下都是0,那么原式的结果就等于给那个数/式子取模后的新式子的结果。

      同余与等式常常转换。

      

  • 相关阅读:
    Unable to open the physical file xxxx. Operating system error 2
    Apache Solr 实现去掉重复的搜索结果
    solr调用lucene底层实现倒排索引源码解析
    lucene-solr本地调试方法
    lucene实战--打分算法没有那么难!
    avalon新一代UI库发布
    让 IE6支持max-height
    迷你MVVM框架 avalonjs 1.3.2发布
    less gradient-vertical 方法的实现
    sass compact方法的实现
  • 原文地址:https://www.cnblogs.com/InductiveSorting-QYF/p/15167925.html
Copyright © 2011-2022 走看看