geometry.md

Вычислительная геометрия

Напомним, что отрезок, для которого указано, какой из его концов считается началом, а какой — концом, называется вектором. Вектор на плоскости можно задать двумя числами — его координатами по горизонтали и вертикали.

Помимо очевидных сложения, вычитания и умножения на константу (скаляр — одно число), у векторов можно ввести и свои особенные операции, которые нам упростят жизнь.

Скалярное произведение (англ. dot product) — произведение длин векторов на косинус угла между ними. Для него справедлива следующая формула:

$$ a \cdot b = x_a x_b + y_a y_b $$

Она доказывается муторно и чисто технически, так что мы это делать не будем.

Геометрически, она равна проекции вектора $b$ на вектор $a$, помноженный на длину $а$:

У него есть полезные свойства:

• Скалярное произведение симметрично ($a \cdot b = b \cdot a$).
• Перпендикулярные вектора должны иметь нулевое скалярное произведение.
• Если угол острый, то скалярное произведение положительное.
• Если угол тупой, то скалярное произведение отрицательное.

Векторное произведение (англ. cross product) — произведение длин векторов на синус угла между ними, причём знак этого синуса зависит от порядка операндов. Оно тоже удобно выражается в координатах:

$$ a \times b = x_a y_b - y_a x_b $$

Геометрически, это ориентированный объем параллелограмма, натянутого на вектора $a$ и $b$:

Его свойства:

• Векторное произведение антисимметрично: $a \times b = - (b \times a)$.
• Коллинеарные вектора должны иметь нулевое векторное произведение.
• Если $b$ «слева» от $a$, то векторное произведение положительное.
• Если $b$ «справа» от $a$, то векторное произведение отрицательное.

Вообще говоря, векторное произведение определяется не так. Оно определено как вектор такой же длины, но перпендикулярный обоим исходным векторам. Это имеет применение в трёхмерной геометрии и физике, но нам об этом думать не надо.

Всякие проверки

Благодаря этим свойствам, почти все проверки в геометрии можно описать через них, а не уравнениями.

Принадлежность точки треугольнику. Пусть у нас есть треугольник $ABC$ (заданный против часовой стрелки) и точка $P$. Тогда она должна лежать слева от всех трёх векторов $AB$, $BC$ и $CA$. Это условие задаст пересечение трёх полуплоскостей, которое и будет нужным треугольником.

$$ \text{P лежит внутри ABC} \iff \begin{cases} (B-A) \times (P-A) \geq 0 \\ (C-B) \times (P-B) \geq 0 \\ (A-C) \times (P-C) \geq 0 \\ \end{cases} $$

Площадь треугольника. Можно пользоваться готовыми формулами, а можно и свойством векторного произведения.

$$ V = \frac{1}{2} (B-A) \times (C-A) $$

Площадь произвольного многоугольника. Если многоугольник задан последовательностью вершин в каком-то порядке, то можно считать так: для каждого ребра добавим его ориентированную площадь от начала координат. Какие-то слагаемые будут положительными (которые на последнем слое, а какие-то — отрицательными).

Забудьте о формуле Герона и всегда считайте площади через векторное произведение.

Кстати, из формулы для площади треугольника следует, что площадь любой фигуры будет либо целым числом, либо рациональным с двойкой в знаменателе. Часто в в задачах входные данные целочисленные, и, чтобы оставаться в целых числах, когда мы считаем какую-нибудь площадь, иногда имеет смысл умножить все входные числа на $2$ (см. «точность»).

Проверка на выпуклость. Можно пройтись по сторонам многоугольника и проверять векторным произведением, что мы поворачиваем всегда в одну сторону, то есть для всех последовательных точек $a$, $b$, $c$ проверить, что $(b-a)\times(c-a) > 0$.

Пересекаются ли отрезки.

Уравнение прямой

Прямую можно задать уравнением вида $Ax + By + C = 0$. Полуплоскость можно задать таким же неравенством.

У прямой есть вектор нормали с координатами $(A, B)$. Он перпендиуклярен прямой, а в случае с полуплоскостью $Ax + By + C \geq 0$ будет указывать в сторону самой полуплоскости.

Чтобы найти расстояние от точки $(x_0, y_0)$ до прямой $Ax + By + C = 0$, можно воспользоваться следующей формулой:

$$ d = \frac{|Ax_0+By_0+C|}{\sqrt{A^2+B^2}} $$

Точка пересечения. По сути, найти точку пересечения двух прямых — это то же самое, что и найти точку, которая удовлетворяет обоим условиям их уравнений:

$$ \begin{cases} A_1 x + B_1 y + C_1 = 0 \\ A_2 x + B_2 y + C_2 = 0 \end{cases} \implies \begin{cases} -x = \frac{B_1 y + C_1}{A_1} \\ -x = \frac{B_2 y + C_2}{A_2} \end{cases} \implies \frac{B_1 y + C_1}{A_1} = \frac{B_2 y + C_2}{A_2} \implies y = - \frac{A_1 C_2 - A_2 C_1}{A_1 B_2 - A_2 B_1} $$

Аналогично, $x = \frac{B_1 C_2 - B_2 C_1}{A_1 B_2 - A_2 B_1}$ (обратите внимание на знаки).

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

Как это кодить в C++

Небольшой ликбез по объектно-ориентированному программированию в C++. Создадим класс, который будет отвечать за все операции с точками. В C++ есть два способа это сделать: через struct и через class. Их основное отличие в том, что по умолчанию в class все поля приватные — к ним нет прямого доступа снаружи. Это нужно для дополнительной защиты, чтобы в крупных промышленных проектах никто случайно ничего не поломал, но на олимпиадах это не очень актуально.

Точка $\simeq$ вектор. Будем считать точка и вектор это один и тот же объект, так как они оба — это просто пара чисел. Будем сопоставлять точке её радиус-вектор — вектор из начала координат, ведущий в эту точку. По принятой в математике и физике нотации, будем обозночать вектора как r. Вы можете обозвать их как point, pt, vec — как угодно.

struct r {
    double x, y;
    r() {}
    r(int _x, int _y) { x = _x, y = _y; }
};

Функция r внутри класса вызывается при инциализации объекта. Её называют конструктором, и её можно указывать разную для разных параметров. Здесь r()вернёт точку с неопределенными (какие оказались в памяти в тот момент) координатами, а r(x, y) вернет точку с координатами $(x, y)$.

Операции над векторами

Давайте напишем функцию, которая принимает вектора и что-то с ними делает. Например, считает длину:

double len(r a) { return sqrt(a.x*a.x + a.y*a.y); }

Операторы

В C++ можно перегружать почти все стандартные операторы, например, +, -, << и т. д.

Переопределим для будущих нужд + и -:

r operator+(r a, r b) { return r(a.x+b.x, a.y+b.y); }
r operator-(r a, r b) { return r(a.x-b.x, a.y-b.y); }

Скалярное произведение:

int operator*(r a, r b) { return a.x*b.x + a.y*b.y; }

Векторное произведение:

int operator^(r a, r b) { return a.x*b.y - b.x*a.y; }

Ввод-вывод

Как вы думаете, как на самом деле работает cin >> x? Это тоже перегрузка оператора — >>. Делается это так:

istream& operator>>(istream &in, r &p) { 
    in >> p.x >> p.y;
    return in;
}

ostream& operator<<(ostream &out, r &p) { 
    out << p.x << " " << p.y << endl;
    return out;            
}

Почему алгебра это плохо

Мы могли не создавать никаких структур и работать с уравнениями, описывающими геометрические объекты. Такой подход будет популярен на олимпиадах по математике, а не по программированию. Когда математик говорит «пересечем две прямые», он представляет громоздкое уравнение, с которым он потом будет работать.

Программист же хочет абстрагироваться и просто написать intersect(a, b), в корректности которого он точно уверен. Программист хочет разбить задачу на много маленьких кусочков и делать по отдельности, а не возиться с формулами.

Приведем несколько примеров конструктивного подхода.

Векторное представление прямой

Прямую можно задать не через уравнение, а через два вектора $a$ и $b$:

$$ Ax + By + C = 0 \rightarrow r = at + b $$

Чтобы это сделать, достаточно выбрать две любые точки на прямой:

// даны A, B, C (A^2 + B^2 != 0)
r a, b;
if (eq(A, 0)) // значит, это горизонтальная прямая
    a = r(0, -C/B), b = r(1, -C/B);
else
    a = r(-C/A, 0), b = (1, -(C+B)/A, 1)

Отражение от прямой

Пусть нам надо отразить точку $(x_0, y_0)$ симметрично относительно заданной прямой $ax+by+c=0$. Чисто в педагогических целях, начнём решать эту задачу как математики, чтобы никогда потом так не делать.

$$ \Pr_a b = \frac{a \cdot b}{|a|} \frac{a}{|a|} = \frac{|a| |b| \cos \alpha}{|a|} \frac{a}{|a|} = |b| \cos \alpha \frac{a}{|a|} $$

Геометрический смысл: длина на единичный вектор направления.

Мы не хотим раскрывать эти формулы покоординатно и предъявлять готовый ответ. Мы знаем, что он получится громоздким. Нам не жалко посчитать всё по частям — здесь нет смысла заниматься оптимизациями. Также мы хотим делать всё по частям, потому что так становится более наглядной логика алгоритма, и, как следствие, его проще дебажить.

// прямая r = at + b, точка c
r pr (r a, r b, r c) {
    c -= b; // пусть c и a выходят из одной точки
    return b + (a*b / len(a) / len(a)) * a;
}

r reflect (r a, r b, r c) {
    return c + 2*(pr(a, b, c)-c);
}

Типичные баги

Точность

Первое правило действительных чисел — не использовать действительные числа

Все переменные типа double хранятся в компьютере неточно (ну а как вы представите ⅓ в двоичной системе счисления?). Поэтому при работе с даблами нужно всегда учитывать эту погрешность. Например, чтобы сравнить два дабла, надо проверить, что они отличаются по модулю меньше, чем на очень маленькое число eps:

const double eps = 1e-8;

bool eq (double a, double b) { return abs(a-b) < eps }

Чтобы так не делать, старайтесь по возможности использовать только инты и абсолютную точность. Иногда есть трюки, позволяющие так делать: например, если в задаче все входные точки целочисленные и нас просят посчитать какую-то площадь, то можно все координаты домножить на два, и тогда ответ тоже будет целым (см. векторное произведение), который только при выводе нужно будет поделить на четыре.

$0 \neq -0$

Действительные числа так хранятся, что $0$ и $-0$ могут быть разными числами. Имейте это ввиду.

Область определения обратных функций

acos, asin и прочие обратные тригонометрические функций требуют, чтобы им на вход подавалось число от -1 до 1. Для безопасности, отмасштабируйте числа, перед тем как брать от них эти функции.