给一个二元一次不定方程
a x + b y = c ax+by=c ax+by=c
- 无整数解:输出 -1
- 有整数解:
- 有正整数解:输出解的个数,输出 x , y x,y x,y 的最小值和最大值
- 无正整数解:输出 x , y x,y x,y 的最小正整数值
思路:
首先是欧几里得算法,用来求 g c d ( a , b ) \mathrm{gcd}(a,b) gcd(a,b),首先假设 a > b a>b a>b,假如 b ∣ a b\mid a b∣a,则 g c d ( a , b ) = b \mathrm{gcd}(a,b)=b gcd(a,b)=b,假如 b ∤ a b\nmid a b∤a,则:
a = k b + c a=kb+c a=kb+c
令 d ∣ a , d ∣ b d\mid a,d\mid b d∣a,d∣b,则显然有 d ∣ c d\mid c d∣c ,令 d = g c d ( a , b ) d=\mathrm{gcd}(a,b) d=gcd(a,b),则 g c d ( a , b ) ∣ c \mathrm{gcd}(a,b)\mid c gcd(a,b)∣c,而 c = a m o d b c=a\bmod b c=amodb,因此 g c d ( a , b ) = g c d ( b , a m o d b ) \mathrm{gcd}(a,b)=\mathrm{gcd}(b,a\bmod b) gcd(a,b)=gcd(b,amodb)。
int gcd(int a,int b){
if(!b)return a;
return gcd(b,a%b);
}
然后是扩展欧几里得,用来求 a x + b y = g c d ( a , b ) ax+by=\mathrm{gcd}(a,b) ax+by=gcd(a,b) 的一组可行解,首先令
a x 1 + b y 1 = g c d ( a , b ) b x 2 + ( a m o d b ) y 2 = g c d ( b , a m o d b ) \begin{aligned} ax_1+by_1&=\mathrm{gcd}(a,b)\\ bx_2+(a\bmod b)y_2&=\mathrm{gcd}(b,a\bmod{b}) \end{aligned} ax1+by1bx2+(amodb)y2=gcd(a,b)=gcd(b,amodb)
根据欧几里得算法,有 g c d ( a , b ) = g c d ( b , a m o d b ) \mathrm{gcd}(a,b)=\mathrm{gcd}(b,a\bmod b) gcd(a,b)=gcd(b,amodb),因此:
a x 1 + b y 1 = b x 2 + ( a m o d b ) y 2 a x 1 + b y 1 = b x 2 + ( a − b ⌊ a b ⌋ ) y 2 a x 1 + b y 1 = a y 2 + b ( x 2 − ⌊ a b ⌋ y 2 ) \begin{aligned} ax_1+by_1&=bx_2+(a\bmod b)y_2\\ ax_1+by_1&=bx_2+(a-b\lfloor\frac{a}{b}\rfloor)y_2\\ ax_1+by_1&=ay_2+b(x_2-\lfloor\frac{a}{b}\rfloor y_2)\\ \end{aligned} ax1+by1ax1+by1ax1+by1=bx2+(amodb)y2=bx2+(a−b⌊ba⌋)y2=ay2+b(x2−⌊ba⌋y2)
可得:
x 1 = y 2 y 1 = x 2 − ⌊ a b ⌋ y 2 \begin{aligned} x_1&=y_2\\ y_1&=x_2-\lfloor\frac{a}{b}\rfloor y_2 \end{aligned} x1y1=y2=x2−⌊ba⌋y2
通过递归可解:
int exgcd(int a,int b,int &x,int &y){
if(!b){x=1,y=0;return a;}
int d=exgcd(b,a%b,x,y);
int t=x;x=y;y=t-(a/b)*y;
return d;
}
由此可以解出 x 0 , y 0 x_0,y_0 x0,y0 使得
a x 0 + b y 0 = g c d ( a , b ) ax_0+by_0=\mathrm{gcd}(a,b) ax0+by0=gcd(a,b)
成立,为了解出 a x + b y = c ax+by=c ax+by=c 的一个特解,对上式做一些变换,得:
a x 0 c g c d ( a , b ) + b y 0 c g c d ( a , b ) = c a\frac{x_0c}{\mathrm{gcd}(a,b)}+b\frac{y_0c}{\mathrm{gcd}(a,b)}=c agcd(a,b)x0c+bgcd(a,b)y0c=c
因此可得原方程得一个特解 x 1 , y 1 x_1,y_1 x1,y1:
x 1 = x 0 c g c d ( a , b ) , y 1 = y 0 c g c d ( a , b ) x_1=\frac{x_0c}{\mathrm{gcd}(a,b)},\;y_1=\frac{y_0c}{\mathrm{gcd}(a,b)} x1=gcd(a,b)x0c,y1=gcd(a,b)y0c
然后构造原方程通解,令 d ∈ Q d\in\mathbb{Q} d∈Q,则
a ( x 1 + d b ) + b ( y 1 − d a ) = c a(x_1+db)+b(y_1-da)=c a(x1+db)+b(y1−da)=c
并且 d a , d b ∈ Z da,db\in\Z da,db∈Z,易知 d = 1 g c d ( a , b ) d=\dfrac{1}{\mathrm{gcd}(a,b)} d=gcd(a,b)1,令
d x = b g c d ( a , b ) , d y = a g c d ( a , b ) d_x=\frac{b}{\mathrm{gcd}(a,b)},\;d_y=\frac{a}{\mathrm{gcd}(a,b)} dx=gcd(a,b)b,dy=gcd(a,b)a
则通解为:
x = x 1 + s d x y = y 1 − s d y \begin{aligned} x&=x_1+sd_x\\ y&=y_1-sd_y \end{aligned} xy=x1+sdx=y1−sdy
其中 s ∈ Z s\in \Z s∈Z。
考虑解出来的 x , y x,y x,y 均为正整数,则
x 1 + s d x > 0 ⟹ s > − x 1 d x y 1 − s d y > 0 ⟹ s < y 1 d y \begin{aligned} x_1+sd_x>0&\implies s>-\frac{x_1}{d_x}\\ y_1-sd_y>0&\implies s<\frac{y_1}{d_y} \end{aligned} x1+sdx>0y1−sdy>0⟹s>−dxx1⟹s<dyy1
因此
L = ⌈ − x 1 + 1 d x ⌉ ≤ s ≤ ⌊ y 1 − 1 d y ⌋ = R L=\lceil\frac{-x_1+1}{d_x}\rceil\le s\le \lfloor\frac{y_1-1}{d_y}\rfloor=R L=⌈dx−x1+1⌉≤s≤⌊dyy1−1⌋=R
假如 L > R L>R L>R,则不存在正整数解, s = L s=L s=L 时得 x min x_{\min} xmin, s = R s=R s=R 时得 x max x_{\max} xmax;
否则,存在正整数解,解的个数是 R − L + 1 R-L+1 R−L+1, s = L s=L s=L 时可以得到 x min , y max x_{\min},y_{\max} xmin,ymax , s = R s=R s=R 时可以得到 x max , y min x_{\max},y_{\min} xmax,ymin
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
typedef long long ll;
int T;
ll a,b,c;
ll exgcd(ll a,ll b,ll &x,ll &y){
if(!b){x=1,y=0;return a;}
ll d=exgcd(b,a%b,x,y);
ll t=x;x=y;y=t-a/b*y;
return d;
}
int main(){
#ifdef WINE
freopen("data.in","r",stdin);
#endif
scanf("%d",&T);
while(T--){
ll x0,y0;
scanf("%lld%lld%lld",&a,&b,&c);
ll d=exgcd(a,b,x0,y0);
if(c%d){
printf("-1\n");
continue;
}
ll x1=x0*c/d,y1=y0*c/d;
ll dx=b/d,dy=a/d;
ll L=(ll)ceil((-x1+1)*1.0/dx);
ll R=(ll)floor((y1-1)*1.0/dy);
ll xmin=x1+L*dx;
ll ymin=y1-R*dy;
if(L>R){
printf("%lld %lld\n",xmin,ymin);
continue;
}
ll xmax=x1+R*dx,ymax=y1-L*dy;
printf("%lld %lld %lld %lld %lld\n",R-L+1,xmin,ymin,xmax,ymax);
}
return 0;
}