浅谈扩欧及exgcd对二元不定方程求解问题

来源:互联网 发布:如何在淘宝上搜索店铺 编辑:程序博客网 时间:2024/06/06 14:24

今天模拟赛第一题是一道对欧几里得和扩欧算法的简单应用,可惜两种方法都不会的我只能用求导和最小矩阵来存不定方程在坐标系上的整数解,满打满算了七十多行代码,其实一个扩欧就能解决的问题,被我想的很复杂。所以,这就是数学结合信息学的恶心之处吧,所以多学学算法对自己还是有帮助的!(博主苦心孤诣地教导各位同学,不要像博主一样碰到数论题就懵逼)

总结一下题意就是求解:ax+by=ca,b,c均为整数)有多少种使得解x,y均为正整数解的情况。首先如果你没有想到扩欧或别的什么算法,直接大暴力是肯定会TLE的,因此我们就要接触一下什么叫做欧几里得和扩欧定理。

欧几里德算法又称辗转相除法,用于计算两个整数a,b的最大公约数。其计算原理依赖于下面的定理:gcd函数就是用来求(a,b)的最大公约数的。
gcd函数的基本性质:
gcd(a,b)=gcd(b,a) 性质一
gcd(a,b)=gcd(-a,b) 性质二
gcd(a,b)=gcd(|a|,|b|) 性质三
gcd(a,b)=gcd(b,a % b) 性质四

对性质的基本证明:
性质一证明:a和b的最大公约数与b和a的最大公约数肯定相等啊
性质二,三证明:若a和b中有一数为负数时,求其gcd时一样求,最后出答案时把“-”丢了就行了,所以负的情况和绝对值情况其实并没多大区别(同一)。
性质四:设d=gcd(a,b),则有d|a,d|b,设r=a%b,因此有b*k+r=a,则r=a-b*k,因此d|r(b*k+a仍为d的倍数),所以则有d|b&&d|r,因此d=gcd(b,a%b),因此gcd(a,b)=gcd(b,a%b)。(这个方法很重要,是求证扩欧的基本应用)

扩展欧几里德算法是用来在已知a, b求解一组x,y使得ax+by = Gcd(a, b) =d(解一定存在,根据数论中的相关定理)。扩展欧几里德常用在求解模线性方程及方程组中。

开始求解:
对于ax+by=c的不定方程,设r=gcd(a,b),
if(c%r!=0)则说明c中不含a,b的最大公约数,因此无整数解
if(c%r==0)将方程右边的c变为c*(r/c)后,变为ax+by=r,此时得到一组整数解:x0,y0(这是转换后方程的一组解,若要求得原解则需把x0,yo都*r/c)
—————————————–手动分割线,滑稽————————————————–
则得到原方程解为x1=x0*(c/r),y1=y0*(c/r);——I
那么x0和y0怎么求解呢?利用我们的性质四不难得出
gcd(a,b)=gcd(b,a%b)
则会有:a*x0+b*y0=b*x1+(a%b)*y1
那么化简得:a*x0+b*y0=b*x1+(a-a/b*b)*y1
继续:a*x0+b*y0=b*x1+a*y1-a/b*b*y1
再继续:a*x0+b*y0=a*y1+b*(x1-a/b*y1)
—————————————–手动分割线,滑稽————————————————–
于是乎我们得到通解为x0=y1,y0=x1-a/b*y1——II
//我们是不是得到了一个递推式了呢?
此处让读者自己思考一下联立I和II,那么就可以得到一个更优化的方程
—————————————–手动分割线,滑稽————————————————–
x=x0+b/r*k,y=y0-a/r*k;
证明如下:对方程两边ax+by=c同除gcd(a,b),则得到Ax+By=C,此时A,B互质(之前的a,b不一定互质的),设r=gcd(a,b)。设此时方程解为x=x0,y=y0(由I得到)代回原方程得Ax0+By0=C;我们为了得到通解,此时设x,y的通解为x=x0+k1,y=y0+k2(k1,k2均为整数),则A(x0+k1)+B(y0+k2)=C,得到Ak1+Bk2=0,则得到Ak1=-Bk2,由于此时A,B互质,所以要使等式成立,k1中应含有因式B,k2中应含有因式A,则设Ak1=-Bk2=A*B*k,则有k1=B*k,k2=-A*k,因此x=x0+B*k,y=y0-A*k,而a/r=A,b/r=B,因此原方程的解为:x=x0+b/r*k,y=y0-a/r*k。好的,大功告成了!我们的解这题的exgcd也引入完毕了(如果读者们觉得累了的话,休息一会儿吧)
—————————————–手动分割线,滑稽————————————————–
那读者们是不是想,只要我们得到第一个状态的x,y值,那后面就好求了呢?那请读者思考一个问题,经过若干次变化后,b是不是会被mod的只剩下0了,那这就是我们所需要的第一组状态,那么在b=0时,gcd(a,b)=a(任何数和0的最大公约数都是它本身),此时x=1,y=0,那么我们就能用递归解决这个问题了!先递归进入下一层,等到到达最后一层即 b=0时就返回x=1,y=0再根据我们的通解得到当层的解不断算出当层的解并返回,最终返回至第一层,得到原解
//好了,现在知识点引入完毕,我们得出通解了(如果对此还有怀疑的同学,可以把我们求出的通解带入最开始的ax+by=c的式子里,我们就发现自己写的是对的了!)

好了,今天讲的这个东西可能比较绕,读者们一定要好好理解,看不懂就多看几遍,这是很重要的,以后肯定会用到的!(其实博主也是翻了百度百科等众多资料后才总结出来的这些精华)。好了,话不多说,现在引入题目

【题目描述】 给出一个二元一次方程 ax+by=c,其中 x、y 是未知数,求它的正整数解的 数量。

【输入格式】 第一行一个整数 T,表示有 T 组数据。接下来 T 行,每行 3 个整数 a、b、c。

【输出格式】 输出 T 行,每行一个数,表示方程解的数量。如果正整数解的数量比 65535 还多,输出“Too much”。

【样例输入】
3
-1 -1 -3
1 1 65536
1 1 65537

【样例输出】
2
65535
Too much

【数据规模与约定】
20%的数据,a=b=1
40%的数据,T≤100,1≤a,b,c≤1000
另20%的数据,a+b=c,1≤a,b,c≤1,000,000
另 20%的数据,1≤a,b,c≤1,000,000
100%的数据,T≤10000,-1,000,000≤a,b,c≤1,000,000

—————————————–手动分割线,滑稽————————————————–
如果读者们已经想清楚了这个GCD的原理的话,那我就上代码了(看不懂一定要反复看,若有疑问或者博主写的有不恰当的地方,欢迎大家指出错误并一起来讨论)
—————————————–手动分割线,滑稽————————————————–
哦,对代码中还有一个小优化,ax+by=c,by=c-ax,y=-a/b*x+c/b,是不是就变成了y=kx+b的形式了?是不是很神奇?那么x,y的解集情况就可以直接在直角坐标系中表现出来了!如果k>0时,y随着x的增大而增大,因此只要求出一组解那么必定会有无穷多组解(x,y都是能增长到正无穷的),如果k<0,b<0时,是没有满足题意的解的,如果k<0,b>0时,我们枚举从x的最小值到x的最大值就行了!

#include <iostream>#include <cstdio>#include <algorithm>#include <cstring>#include <cmath>using namespace std;#define lol long longlol gcd(lol x,lol y){    return y?gcd(y,x%y):x;}void exgcd(lol a,lol b,lol &x,lol &y){    if(b==0)    {        x=1;        y=0;        return;    }    exgcd(b,a%b,y,x);    y=y-a/b*x;}lol a,b,c,d,T;int main(){    scanf("%lld",&T);    while(T--)    {        scanf("%lld%lld%lld",&a,&b,&c);        if(a==0||b==0)        {            if(a==0)swap(a,b);            if(a==0)            {                if(!c) printf("ZenMeZheMeDuo\n");                else printf("0\n");            }            else            {                lol d=c/a;                if(d>0&&a*d==c) printf("ZenMeZheMeDuo\n");                else printf("0\n");            }            continue;        }        //特判处理完毕        lol r=gcd(a,b);        lol k1=b/r,k2=a/r,x,y,z,ans;        if(c%r!=0)        {            printf("0\n");            continue;        }        exgcd(a,b,x,y);        x*=c/r;        y*=c/r;        if(k1<0)        {            k1=-k1;            k2=-k2;        }        if(a*b<0) printf("ZenMeZheMeDuo\n");        else        {            if(c*b<0) printf("0\n");            else            //解出的x0,y0的第一组解可能并不是我们所想要的,比如x0y0有可能为负数,此时我们需要将他们变为正数才方便求解            {                if(x<0)                {                    z=-x/k1;                    x+=z*k1;                    y-=z*k2;                    while(x<=0)                    {                        x+=k1;                        y-=k2;                    }                }                else                {                    z=x/k1;                    x-=z*k1;                    y+=z*k2;                    while(x<=0)                    {                        x+=k1;                        y-=k2;                    }                }                if(y<=0) ans=0;                else ans=y/k2+bool(y%k2);                if(ans>65535) printf("ZenMeZheMeDuo\n");                else printf("%lld\n",ans);            }        }    }    return 0;}

写到这里,码了两天终于写完了,跑了24ms,感动ing,希望自己能更多接触一些数论方面的东西,否则到时优化算法时会很吃亏的。大家一起勉励,一起加油,一起讨论!
Thanks for your attention.