Пересечение двух лучей

Спасибо @MBo за контрпример. Действительно, заплатка в виде a == c, призванная поддерживать случаи по типу:

0 0
0 1
0 0
0 -1

меня напрягала. На самом деле, достаточно проверить, лежит ли начало у каждого луча на другом луче при совпадении прямых, т.е.:

// Предполагаем, что точка точно лежит на прямой
bool is_point_on_the_ray(Point a, Point b, Point p) {
    return scalar_prod(Vector(a, b), Vector(a, p)) >= 0;
}

int main() {
    // ...
    if (line_intersect_status == 2) {
        are_rays_intersect = is_point_on_the_ray(c, d, a) || is_point_on_the_ray(a, b, c);
    } else if (line_intersect_status == 1) {
        are_rays_intersect = is_point_on_the_ray(a, b, p) && is_point_on_the_ray(c, d, p);
    } else are_rays_intersect = false;
    
    // ...
}

Но опять же, отмечать это как ответ я не буду, возможно, мне в какой-то степени повезло, что точности long double хватило. В комментариях написали, что есть более элегантный способ без использования деления, поэтому я отдам все лавры ему, если таковой ответ здесь появится.

Даны три точки: p₁ = (x₁, y₁), p₂ = (x₂, y₂) и p = (x, y). Тогда определитель ниже равен удвоенной площади треугольника Δp₁p₂p (общая формула площади треугольника в декартовых координатах на плоскости). Площадь со знаком, она положительна если точки расположены против часовой стрелки, отрицательна если точки - по часовой стрелке, равна нулю, если точки на одной прямой.

Зафиксируем p₁, p₂. Будем менять только p. Определить положителен если p "слева" от прямой p₁p₂, отрицателен если "справа", равен нулю если точка на прямой. Значение определителя - линейная функция от координат x и y:

d(x, y) = (y₁ - y₂)x + (x₂ - x₁)y + (x₁y₂ - x₂y₁).

Луч задаётся парой точек: p - начало луча, q - точка на луче. Любую точку r на луче можно представить в виде r(t) = (1 - t)p + tq, где t - неотрицательное число. В частности r(0) = p, r(1) = q. В координатах будет:

x_t = (1 - t)p_x + tq_x, y_t = (1 - t)p_y + tq_y, где t ≥ 0.

Обозначим f(t) = d(x_t, y_t). d, x_t, y_t - линейные функции, следовательно f тоже линейная.

Луч r(t) пересекается с прямой p₁p₂ тогда и только тогда, когда существует t ≥ 0 такое что f(t) = 0.

f - линейная. Чтобы определить что она имеет корень, достаточно проанализировать значения f(0) и f(1). Например, если 0 < f(0) ≤ f(1), корня точно нет. Тоже верно для случая f(1) ≤ f(0) < 0. Если f(0) = f(1) = 0, луч целиком лежит на прямой. Иначе луч пересекается с прямой в одной точке.

Вся эта теория используется в функции int intersect(const LineEq &eq, const Ray &ray), которая определяет как размещены на плоскости прямая и луч. Они могут не пересекаться совсем, пересекаться в одной точке, луч может лежать на прямой. Функция f явно не строится, нам нужны только значения f(0) = d(p_x, p_y) и f(1) = d(q_x, q_y).

Функция bool intersect(const Ray &r1, const Ray &r2) определяет пересечение двух лучей:

• если первый луч не пересекает вторую прямую (продолжение второго луча), то и лучи не пересекаются;
• если второй луч не пересекает первую прямую, лучи не пересекаются;
• если оба пересечения луч-прямая происходят в одной точке, то и сами лучи пересекаются в этой точке;
• если оба пересечения луч-прямая происходят во многих точках, лучи лежат на общей прямой, этот случай требует дополнительной работы.

Пусть луч задан началом p = (p_x, p_y) и точкой на луче q = (q_x, q_y). На прямой pq лежит точка r = (x, y). Тогда r лежит на луче, тогда и только тогда когда (r - p)(q - p) ≥ 0. В координатах:

(x - p_x)(q_x - p_x) + (y - p_y)(q_y - p_y) = (q_x - p_x)x + (q_y - p_y)y + (p_x(p_x - q_x) + p_y(p_y - q_y)).

Последнее выражение линейная функция g(x, y). Подставим в неё координаты точек из второго луча и проверим принимает ли она неотрицательные значения. Логика очень похожа на пересечение луча и прямой, только проще.

Это вся теория:

#include <iostream>

struct Point { int x, y; };
struct Ray { Point p, q; };

// |x1 y1 1|
// |x2 y2 1| = (y1 - y2)x + (x2 - x1)y + (x1y2 - x2y1)
// |x  y  1|
class LineEq {
public:
    LineEq(const Point &p1, const Point &p2)
    // 2M             2M              2MM
    : a(p1.y - p2.y), b(p2.x - p1.x), c(p1.x * p2.y - p2.x * p1.y) {}
    //                                            6MM
    int operator()(const Point &p) const { return a * p.x + b * p.y + c; }
private:
    int a, b, c;
};

// 0 - no intersection, 1 - one point intersection, 2 - ray lays on line
int intersect(const LineEq &eq, const Ray &ray) {
    const int f0 = eq(ray.p);
    const int f1 = eq(ray.q);
    if (f1 <= f0 && f0 < 0) { return 0; }
    if (0 < f0 && f0 <= f1) { return 0; }
    return (f0 != f1) ? 1 : 2;
}

// (r - p)(q - p) =
// = (x - px)(qx - px) + (y - py)(qy - py) =
// = (qx - px)x + (qy - py)y + (px(px - qx) + py(py - qy))
class RayEq {
public:
    RayEq(const Point &p, const Point &q)
    // 2M             2M              4MM
    : a(q.x - p.x), b(q.y - p.y), c(p.x * -a + p.y * -b) {}
    //                                            8MM
    int operator()(const Point &r) const { return a * r.x + b * r.y + c; }
private:
    int a, b, c;
};

bool intersect(const Ray &r1, const Ray &r2) {
    const int i12 = intersect(LineEq(r1.p, r1.q), r2);
    if (i12 == 0) {
        return false;
    }
    const int i21 = intersect(LineEq(r2.p, r2.q), r1);
    if (i21 == 0) {
        return false;
    }
    if (i12 == 1) {
        return true;
    }

    const RayEq eq(r1.p, r1.q);
    const int f0 = eq(r2.p);
    if (0 <= f0) {
        return true;
    }
    const int f1 = eq(r2.q);
    return f0 < f1;
}

int main() {
    Ray r1, r2;
    if (std::cin
        >> r1.p.x >> r1.p.y
        >> r1.q.x >> r1.q.y
        >> r2.p.x >> r2.p.y
        >> r2.q.x >> r2.q.y
    ) {
        std::cout << ((intersect(r1, r2)) ? "YES" : "NO") << '\n';
    }
}

Решение использует только целые числа и операции +, -, * и сравнения. Если модуль координат точек не превосходит M, то максимальное значение в коде не превосходит 8M². Соответствующие комментарии есть в коде. Для типа int требуется 8M² < 2³¹. То есть, переполнения не будет пока модуль координат меньше 2¹⁴.