使用扩展欧几里得算法查找线段与二维网格上的点的交点数答案

【问题标题】：Using extended euclidean algorithm to find number of intersections of a line segment with points on a 2D grid使用扩展欧几里得算法查找线段与二维网格上的点的交点数
【发布时间】：2014-09-02 21:08:40
【问题描述】：

在由网格点 (M*x, M*y) 构造的网格中，给定点 A(x1,y1) 和点 B(x2,y2)，其中所有变量都是整数。我需要检查从 A 点到 B 点的线段上有多少网格点。我知道可以通过某种方式使用扩展欧几里得算法来完成，但我不知道该怎么做。非常感谢您的帮助。

【问题讨论】：

标签： algorithm discrete-mathematics greatest-common-divisor number-theory

【解决方案1】：

换个说法，你要确定t满足多少个数字

(1) M divides (1-t) x1 + t x2
(2) M divides (1-t) y1 + t y2
(3) 0 <= t <= 1.

让我们关注(1)。我们引入一个整数变量q来表示整除性约束并求解t：

exists integer q,  M q = (1-t) x1 + t x2

exists integer q,  M q - x1 = (x2 - x1) t.

如果x1 不等于x2，则这给出了t in {a + b q | q integer} 形式的周期性可能性集，其中a 和b 是分数。否则，所有t 或没有t 都是解决方案。

扩展的欧几里得算法对于与(1) 和(2) 产生的解集相交非常有用。假设我们要计算交集

{a + b q | q integer} intersect {c + d s | s integer}.

通过用两种不同的方式表达一个通用的公共元素，我们得到了linear Diophantine equation：

a + b q = c + d s,

其中a, b, c, d 是常量，q, s 是整数。让我们清除分母并将项收集到一个等式中

A q + B s = C,

其中A, B, C 是整数。这个方程有解当且仅当A 和B 的最大公约数g 也整除C。使用扩展欧几里得算法计算整数系数u, v 使得

A u + B v = g.

那么我们就有办法了

q = u (C/g) + k (B/g)
s = v (C/g) - k (A/g)

对于每个整数k。

最后，我们必须考虑约束(3)。这应该归结为一些算术和一层除法，但我宁愿不计算细节（到目前为止我所写的特殊情况已经花费了你很多时间）。

【讨论】：

【解决方案2】：

让我们dX = Abs(x2-x1) 和dY = Abs(y2 - y1)
那么段上的格点数为

P = GCD(dX, dY) + 1

（包括起点和终点）
其中 GCD 是最大公约数（通过通常的（未扩展的）欧几里得算法）

（查看最后一个属性here）

【讨论】：

【解决方案3】：

按照 David Eisenstat 先生的指示，我设法用 c++ 编写了一个计算答案的程序。

#include <iostream>
#include <math.h>

using namespace std;

int gcd (int A, int B, int &u, int &v)
{
int Ad = 1;
int Bd = 1;
if (A < 0) { Ad = -1; A = -A; }
if (B < 0) { Bd = -1; B = -B; }

int x = 1, y = 0;
u = 0, v = 1;

while (A != 0)
{
    int q = B/A;
    int r = B%A;
    int m = u-x*q;
    int n = v-y*q;
    B = A;
    A = r;
    u = x;
    v = y;
    x = m;
    y = n;
}

u *= Ad;
v *= Bd;

return B;
}

int main(int argc, const char * argv[])
{
int t;
scanf("%d", &t);

for (int i = 0; i<t; i++) {

    int x1, x2, y1, y2, M;

    scanf("%d %d %d %d %d", &M, &x1, &y1, &x2, &y2);

    if ( x1 == x2 ) // vertical line
    {
        if (x1%M != 0) { // in between the horizontal lines
            cout << 0 << endl;
        } else
        {
            int r = abs((y2-y1)/M); // number of points
            if (y2%M == 0 || y1%M == 0) // +1 if the line starts or ends on the point
                r++;

            cout << r << endl;
        }

    } else {

        if (x2 < x1)
        {
            int c;
            c = x1;
            x1 = x2;
            x2 = c;
        }

        int A, B, C;

        C = x1*y2-y1*x2;
        A = M*(y2-y1);
        B = -M*(x2-x1);

        int u, v;
        int g = gcd(A, B, u, v);

        //cout << "A: " << A << " B: " << B << " C: " << C << endl;
        //cout << A << "*" << u <<"+"<< B <<"*"<<v<<"="<<g<<endl;

        double a = -x1/(double)(x2-x1);
        double b = M/(double)(x2-x1);

        double Z = (-a-C*b/g*u)*g/(B*b);
        double Y = (1-a-C*b/g*u)*g/(B*b);

        cout << floor(Z) - ceil(Y) + 1 << endl;
    }
}
return 0;
}

感谢您的帮助！请检查是否考虑了所有特殊情况。

【讨论】：