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

原文传送：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中的一个正确值就好。

　　　　故对于

　　启发思路：

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

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

　　杂想：扩展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*log₂10¹⁸ =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，那么原式的结果就等于给那个数/式子取模后的新式子的结果。

　　同余与等式常常转换。