题解 P5656 【【模板】二元一次不定方程(exgcd)】

· · 题解

2020/1/25 UPD:更正了证明过程中符号的问题。(By weapon_zipper)

2019/12/11 UPD:更正了在题解页面Latex代码漏出的问题。

该题解主要由以下内容组成:

扩展欧几里得算法(exgcd)

扩展欧几里得算法:用于求解形如ax+by=gcd(a,b)的不定方程特解。

求解方法&证明: 当b=0时,可以看出gcd(a,b)=a,而

\left\{\begin{aligned}x & = 1\\y & = 0\\\end{aligned}\right.

是方程的一组特解。

当b\neq0时,递归求解exgcd(b,a%b,x,y),设

\left\{\begin{aligned}a' & = b\\b' & = a\%b\\\end{aligned}\right.

可以求得a'x'+b'y'=gcd(a,b)的一组特解,即x',y'。

因为a'x'+b'y'等价于bx'+(a%b)y',再根据模运算的定义可得bx'+(a-a/b*b)y'

bx'+ay'-(a/b*b)y'

ay'+b(x'-a/b*y')

此时就得到了

\left\{\begin{aligned}x & = y'\\y & = x'-a/b*y'\\\end{aligned}\right.

代码:


void exgcd(int a,int b,int &x,int &y)
{
    if(!b)
    {
        x=1;
        y=0;
        return;
    }
    int p;
    exgcd(b,a%b,x,y);
    p=x;
    x=y;
    y=p-(a/b)*y;
    return;
}

求二元一次方程/线性同余方程特解

二元一次方程: 对于形如ax+by=c的二元一次方程,根据扩展欧几里得算法的定义可得,当且仅当gcd(a,b)|c(|为整除)时,存在整数解。

设g=gcd(a,b),a'=a/g,b'=b/g,c'=c/g,则ax+by=c\Leftrightarrowa'x+b'y=c',

此时gcd(a',b')=1,可以利用exgcd求出a'x'+b'y'=1的一组特解,继而得出

\left\{\begin{aligned}x_0 & = c'x'\\y_0 & = c'y'\\\end{aligned}\right.

我们便求得了ax+by=c的一组特解。

线性同余方程: 对于形如ax\equivc(mod b)的线性同余方程,根据模运算的定义,在方程左侧添加一个by不会对结果造成影响,其实质就等价于ax+by=c的不定方程,利用上述方法求解便可。

求二元一次方程/线性同余方程通解

二元一次方程: 对于形如ax+by=c的二元一次方程,整数解通解为

\left\{\begin{aligned}x_t & = x_0+b't\\y_t & = y_0-a't\\\end{aligned}\right.

(t\in\mathbb{Z})

对其进行证明:ax_{t}+by_{t}

a(x_{0}+b't)+b(y_{0}-a't)

ax_{0}+ab't+by_{0}+a'bt

ax_{0}+by_{0}+at*

\displaystyle\frac{b}{g}$-bt$* \displaystyle\frac{a}{g}

ax_{0}+by_{0}+\displaystyle\frac{abt}{g}-\displaystyle\frac{abt}{g}

ax_{0}+by_{0}

因为lcm(a,b)=\displaystyle\frac{ab}{g},x,y的系数分别为a,b,所以x,y的最小增减幅度分别为\displaystyle\frac{b}{g}\displaystyle\frac{a}{g}即b',a'。

线性同余方程: 同理于二元一次方程。

暴力解决本题

观察本题,发现对于给出的条件有三种情况:无整数解、无正整数解但有整数解、有正整数解。

无整数解: 要求输出-1。根据扩展欧几里得算法的定义可得,当c%gcd(a,b)\ne0时方程无整数解,输出-1即可。

判断有无正整数解: 根据二元一次方程的通解可得,x,y的增减性相反,所以若存在正整数解,则当x取得最小正整数解时,y取得最大正整数解,若此时y\le0,则说明该方程无正整数解。

无正整数解但有整数解: 要求输出x,y的最小正整数解。通过整除运算计算出x取得最小正整数解时的k值,继而求出x的最小正整数解,同理可以求出y的最小正整数解。

有整数解: 此时要求最复杂,要求输出正整数解的数量、x,y的最小正整数解、x,y的最大正整数解。因为当x取得最小正整数解时,y取得最大正整数解,可以利用上述方法求出x取得最小正整数解时的k值,x的最小正整数解与y的最大正整数解,该方法同样可以用于求解x的最大正整数解与y的最小正整数解。正整数解的数量=\displaystyle\frac{x_{max}-x_{min}}{b'}+1或\displaystyle\frac{y_{max}-y_{min}}{a'}+1。

代码:

#include<cstdio>
using namespace std;
int read()
{
    int res=0;
    char c;
    c=getchar();
    while(c<'0'||c>'9')
    c=getchar();
    while(c>='0'&&c<='9')
    {
        res=res*10+c-48;
        c=getchar();
    }
    return res;
}//因为题目描述“读入输出量较大,注意使用较快的读入输出方式”,所以使用快读,实际上不用也可以通过本题。 
long long gcd(long long n,long long m)
{
    return (n%m==0)?m:gcd(m,n%m);
}//gcd,不做说明。 
void exgcd(long long a,long long b,long long &x,long long &y)
{
    if(!b)
    {
        x=1;
        y=0;
        return;
    }
    long long p;
    exgcd(b,a%b,x,y);
    p=x;
    x=y;
    y=p-(a/b)*y;
    return;
}//exgcd,本题核心代码,在上述内容中已做说明 。 
int t;
int main()
{
    t=read();
    while(t--)
    {
        long long x=0,y=0,a,b,c,g,xin,yin,xax,yax,npa=0,k;//t,x,y,a,b,c,g同先前描述,xin,yin为x,y的最小正整数解,xax,yax为x,y的最大正整数解,npa为正整数解数。 
        a=read();
        b=read();
        c=read();//读入部分。 
        g=gcd(a,b);
        if(c%g!=0)
        printf("-1\n");//判断有无整数解。 
        else
        {
            a/=g;
            b/=g;
            c/=g;
            exgcd(a,b,x,y);
            x*=c;
            y*=c;//求出一组特解。 
            for(int i=1;i<=2;i++)//因为在求解y取得最小正整数解时的k值可能使x的值小于等于0,所以要再求解一遍以求出x的最小正整数解与y的最大正整数解。 
            if(x<=0)
            {
                k=-(x/b)+1;
                x+=k*b;
                y-=k*a;
                xin=x;
                yax=y;
            }//求出x取得最小正整数解时的k值,x的最小正整数解与y的最大正整数解。 
            else if(y<=0)
            {
                k=-(y/a)+1;
                x-=k*b;
                y+=k*a;
                yin=y;
                xax=x;
            }//求出y取得最小正整数解时的k值,x的最大正整数解与y的最小正整数解。 
            if(x>0&&y>0)//判断方程有无正整数解。 
            {
                if(x%b!=0)
                {
                    xin=x%b;
                    yax=y+a*(x/b);
                }
                else
                {
                    xin=x%b+b;
                    yax=y+a*(x/b-1);
                }
                if(y%a!=0)
                {
                    yin=y%a;
                    xax=x+b*(y/a);
                }
                else
                {
                    yin=y%a+a;
                    xax=x+b*(y/a-1);
                }//用于求解x,y均大于等于0时x,y的最小最大正整数解。 
                npa=(xax-xin)/b+1;//求出正整数解的数量。 
            }
            if(!npa)
            printf("%lld %lld\n",xin,yin);
            else printf("%lld %lld %lld %lld %lld\n",npa,xin,yin,xax,yax);//输出部分。 
        }
    }
    return 0;
}

对暴力进行优化

由于以上代码不优化显得又长又慢常数巨大而且显得丑,所以尝试对暴力进行优化。

运用大眼观察法,可以发现我们能用模运算直接求出x,y的最小正整数解而不用求k值,对于x,y的最大正整数解则可以直接代入二元一次方程求解,即

\left\{\begin{aligned}x_{max}=\displaystyle\frac{c-x_{min}}{a}\\y_{max}=\displaystyle\frac{c-y_{min}}{b}\\\end{aligned}\right.

代码:

#include<cstdio>
using namespace std;
int read()
{
    int res=0;
    char c;
    c=getchar();
    while(c<'0'||c>'9')
    c=getchar();
    while(c>='0'&&c<='9')
    {
        res=res*10+c-48;
        c=getchar();
    }
    return res;
}
long long gcd(long long n,long long m)
{
    return (n%m==0)?m:gcd(m,n%m);
}
void exgcd(long long a,long long b,long long &x,long long &y)
{
    if(!b)
    {
        x=1;
        y=0;
        return;
    }
    long long p;
    exgcd(b,a%b,x,y);
    p=x;
    x=y;
    y=p-(a/b)*y;
    return;
}
int t;
int main()
{
    t=read();
    while(t--)
    {
        long long x=0,y=0,a,b,c,g,xin,yin,xax,yax,npa=0,k;
        a=read();
        b=read();
        c=read();
        g=gcd(a,b);
        if(c%g!=0)
        printf("-1\n");
        else
        {
            a/=g;
            b/=g;
            c/=g;
            exgcd(a,b,x,y);
            x*=c;
            y*=c;
            xin=x>0&&x%b!=0?x%b:x%b+b;//若x>0且x%b!=0,x的最小正整数解就等于x%b,否则则需要在x%b的基础上加b。 
            yax=(c-xin*a)/b;//求出对应的y的最大正整数解。 
            yin=y>0&&y%a!=0?y%a:y%a+a;
            xax=(c-yin*b)/a;//同理于x。 
            if(xax>0)//判断有无正整数解。 
            npa=(xax-xin)/b+1;//求出正整数解的数量。 
            if(!npa)
            printf("%lld %lld\n",xin,yin);
            else printf("%lld %lld %lld %lld %lld\n",npa,xin,yin,xax,yax);
        }
    }
    return 0;
}

再开个O2,交一发上去,发现跑得还挺快,能拿rank6(主要是因为人少)。 最后一点:记开long long

为什么我会这么不熟......