UVa 11768 Lattice Point or Not (扩展欧几里得)

UVa 11768 Lattice Point or Not

---

题目大意:

给两个点A(x1,y1)和B(x2,y2).其中x1,y1,x2,y2皆为0.1的整数倍,且绝对值不超过200000.统计线段AB经过的整点数.

题目分析:

(看来扩欧我还是不会欸,orz……)

求直线上的整点要用到扩展欧几里得解线性方程.
先利用A,B两点表示出形如$ax+by=c$形式的直线方程,得到
(注意A,B横纵坐标皆为0.1的整数倍,所以事先要将$x1,y1,x2,y2$扩大10倍).

$(y2-y1)x+(x1-x2)y=y2x1-y1x2$
$a=y2-y1,b=x1-x2,c=y2x1-y1x2$
(a,b是坐标相减,c是坐标相乘,所以要将a,b再扩大10倍)

扩展欧几里得可以求解直线$ax+by=gcd(a,b)$的一组整数解$(x0,y0)$.
首先就可以判断,若c不为g倍数,直线无法过任意整点.
若c为g倍数,则有直线$ax+by=c$的一组整数解为$(x0*c/g,y0*c/g)$.

既然已经有了一组整数解,现在要做的就是求出每隔多少,就会又有一组整数解.
设一组整数解为$(x',y')$,则有

$ax+by=ax'+bx'$
$a(x-x')=b(y'-y)$

将$a,b$同除$gcd(a,b)$得$a',b'$,且$a,b$互素,则有

$x-x'=kb',y'-y=ka'$
$x'=x-kb',y'=y+ka'$

此时所求的$a'$与$b'$即为整数解出现的周期,只需要用其中一个即可求解.

以$b'$为例,已知$b'$后,先将$x$值用$b'$改变至[x1,x2]内,再求从$x$到$x2$的解的个数.
(因为$x1,x2$可能不为合法解,不能直接用$(x2-x1)/b+1$求解).

代码:

#include<cmath>
#include<cstdio>
#include<iostream>
#include<algorithm>

using namespace std;

typedef long long ll;

void exgcd(ll a,ll b,ll& g,ll& x,ll& y)
{
    if(!b) g=a,x=1,y=0;
    else exgcd(b,a%b,g,y,x),y-=x*(a/b);
}

double X1,Y1,X2,Y2;

ll solve()
{
    ll x1=(X1+0.05)*10,y1=(Y1+0.05)*10,x2=(X2+0.05)*10,y2=(Y2+0.05)*10;//对拍的时候发现不加0.05貌似会有误差 
    if(x1==x2) {//平行y轴 
        if(x1%10) return 0;
        if(Y2<Y1) swap(Y1,Y2);
        return floor(Y2)-ceil(Y1)+1;
    }
    if(y1==y2) {//平行x轴 
        if(y1%10) return 0;
        if(X2<X1) swap(X1,X2);
        return floor(X2)-ceil(X1)+1;
    }
    ll a=(y2-y1)*10,b=(x1-x2)*10,c=y2*x1-y1*x2,g,x,y;//ax+by=c,注意a,b需要扩大10倍 
    exgcd(a,b,g,x,y);//找到在ax+by=gcd(a,b)条件下的一个合法解(x,y) 
    if(c%g) return 0;//c不为gcd(a,b)倍数,无解 
    x*=c/g;b=abs(b/g);//找到在ax+by=c条件下的合法解的横坐标x 
    if(X1>X2) swap(X1,X2);
    x1=ceil(X1);//小数向上取整 
    x2=floor(X2);//大数向下取整 
    x-=(x-x1)/b*b;//将x改变值至[x1,x2] 
    if(x<x1) x+=b;
    if(x2<x) return 0;
    return (x2-x)/b+1;
}

int main()
{
    int T;
    scanf("%d",&T);
    while(T--) {
        scanf("%lf%lf%lf%lf",&X1,&Y1,&X2,&Y2);
        printf("%lld\n",solve());
    }
    return 0;
}