Вычислительная геометрия: однородные координаты — точки, прямые, пересечения

• Геометрия
• Однородные координаты
• Python
• C/C++
• Rust

Зачем нужны однородные координаты

Точка плоскости обычно задаётся парой $(x, y)$, а прямая — уравнением $ax + by + c = 0$, где $a$, $b$, $c$ — числа и $(a, b) \neq (0, 0)$. Вычисления с такими уравнениями требуют разбора случаев. Пусть даны прямые $y = k_1 x + b_1$ и $y = k_2 x + b_2$. Их общая точка удовлетворяет равенству $k_1 x + b_1 = k_2 x + b_2$, откуда $x = (b_2 - b_1)/(k_1 - k_2)$. Если $k_1 = k_2$, делить нельзя: прямые параллельны или совпадают. Например, для $y = 2x + 1$ и $y = 2x - 3$ знаменатель равен нулю. Кроме того, вертикальная прямая $x = 1$ не записывается в виде $y = kx + b$ ни при каких $k$ и $b$.

Однородные координаты устраняют такие исключения. К паре $(x, y)$ добавляется третья координата $W$, и точка записывается тройкой $(X, Y, W)$. Точки с $W = 0$ называются точками на бесконечности. Тройка $(X, Y, 0)$ задаёт направление $(X, Y)$, и параллельные прямые имеют общую точку на бесконечности. Например, все прямые $2x - y + c = 0$ при любом $c$ проходят через точку $(1, 2, 0)$: подставляем $X = 1$, $Y = 2$, $W = 0$ и получаем $2 \cdot 1 - 1 \cdot 2 + c \cdot 0 = 0$. Благодаря этому любые прямые, в том числе вертикальные и параллельные координатным осям, записываются вектором и обрабатываются одними и теми же формулами.

Точка в однородных координатах

Однородные координаты точки — ненулевая тройка чисел $(X, Y, W)$. Тройки $(X, Y, W)$ и $\lambda (X, Y, W)$ при $\lambda \neq 0$ задают одну и ту же точку; её обозначаем $(X : Y : W)$. Точка плоскости — это, таким образом, класс эквивалентности пропорциональных троек.

Если $W \neq 0$, то точка имеет декартовы координаты $x = X/W$ и $y = Y/W$. Они не зависят от выбора тройки: для $\lambda (X, Y, W)$ получаем $\lambda X/(\lambda W) = X/W$ и $\lambda Y/(\lambda W) = Y/W$. Обратно, точка $(x, y)$ записывается тройкой $(x, y, 1)$, то есть $(x : y : 1)$.

Примеры. Тройки $(2, 4, 2)$, $(1, 2, 1)$ и $(-3, -6, -3)$ задают одну точку $(1, 2)$: в каждом случае $X/W = 1$ и $Y/W = 2$, а $(-3, -6, -3) = -3 \cdot (1, 2, 1)$. Тройка $(3, 6, 0)$ задаёт точку на бесконечности $(3 : 6 : 0) = (1 : 2 : 0)$, то есть направление $(1, 2)$.

Прямая как вектор

Прямую $ax + by + c = 0$ при $(a, b) \neq (0, 0)$ задаём вектором $l = (a, b, c)$. Множество решений не меняется, если все коэффициенты умножить на $\lambda \neq 0$, поэтому вектор $l$ определён с точностью до ненулевого множителя. Скалярное произведение троек определяем формулой $u \cdot v = u_1 v_1 + u_2 v_2 + u_3 v_3$.

Утверждение. Точка $P = (X, Y, W)$ лежит на прямой $l$ тогда и только тогда, когда $P \cdot l = 0$.

Доказательство. Пусть $W \neq 0$, $x = X/W$, $y = Y/W$. Тогда $P \cdot l = aX + bY + cW = W(ax + by + c)$, и так как $W \neq 0$, равенство $P \cdot l = 0$ равносильно $ax + by + c = 0$. Пусть $W = 0$. Тогда $P \cdot l = aX + bY$. В плоскости все векторы, ортогональные ненулевому вектору $(a, b)$, пропорциональны $(-b, a)$, а $(-b, a)$ — направление прямой $l$: вдоль него выражение $ax + by$ не меняется, так как $a \cdot (-b) + b \cdot a = 0$. Поэтому равенство $P \cdot l = 0$ означает, что $(X, Y)$ — направление прямой $l$; именно так точка на бесконечности и должна лежать на прямой.

Пример. Прямая $3x - 2y + 6 = 0$ имеет вектор $l = (3, -2, 6)$. Точка $(0, 3)$ на ней лежит: $(0, 3, 1) \cdot (3, -2, 6) = -6 + 6 = 0$. Точка $(2, 0)$ не лежит: $(2, 0, 1) \cdot (3, -2, 6) = 6 + 6 = 12 \neq 0$. Точка на бесконечности $(2 : 3 : 0)$ лежит на прямой, так как $(2, 3)$ — её направление, и $(2, 3, 0) \cdot (3, -2, 6) = 6 - 6 = 0$.

Прямая через две точки и точка пересечения

Для троек $u = (u_1, u_2, u_3)$ и $v = (v_1, v_2, v_3)$ векторное произведение определяем формулой $$u \times v = (u_2 v_3 - u_3 v_2,\; u_3 v_1 - u_1 v_3,\; u_1 v_2 - u_2 v_1).$$ Раскрыв скобки, получаем тождества $(u \times v) \cdot u = 0$ и $(u \times v) \cdot v = 0$: в выражении $(u \times v) \cdot u$ шесть слагаемых попарно сокращаются, для $v$ рассуждение то же.

Факт. Если $u \neq 0$ и $u \times v = 0$, то $v = \mu u$ для некоторого числа $\mu$.

Доказательство. Выберем $k$ с $u_k \neq 0$. Каждая компонента $u \times v$ равна $\pm(u_i v_j - u_j v_i)$ для некоторой пары $i \neq j$, поэтому из $u \times v = 0$ следует $u_k v_j = u_j v_k$ при всех $j$. Значит, $v_j = (v_k / u_k) \, u_j$ при всех $j$, и можно взять $\mu = v_k / u_k$.

Прямая через две точки. Пусть $P$ и $Q$ — различные точки, из которых хотя бы одна конечна, и $l = P \times Q$. Тогда $l \neq 0$: иначе по факту $Q = \mu P$, и так как $Q \neq 0$, то $\mu \neq 0$, то есть $P$ и $Q$ совпадают как точки. Тождества выше дают $(P \times Q) \cdot P = 0$ и $(P \times Q) \cdot Q = 0$, и по утверждению раздела «Прямая как вектор» обе точки лежат на $l$. Пусть, например, $P$ конечна ($p_3 \neq 0$). Тогда $(l_1, l_2) \neq (0, 0)$: иначе $l = (0, 0, l_3)$, и равенство $P \cdot l = p_3 l_3 = 0$ даёт $l_3 = 0$, то есть $l = 0$. Значит, $l$ — прямая в смысле определения, $(a, b) \neq (0, 0)$. Если же обе точки бесконечно удалённые, то $P \times Q = (0, 0, l_3)$ с $l_3 \neq 0$: это бесконечно удалённая прямая, которая определением не допускается. Поэтому прямая через две точки строится, когда хотя бы одна из точек конечна.

Пример. Пусть $P = (0, 0, 1)$ и $Q = (2, 1, 1)$. Тогда $$l = P \times Q = (0 \cdot 1 - 1 \cdot 1,\; 1 \cdot 2 - 0 \cdot 1,\; 0 \cdot 1 - 0 \cdot 2) = (-1, 2, 0),$$ то есть прямая $-x + 2y = 0$. Проверка: $(0, 0, 1) \cdot (-1, 2, 0) = 0$ и $(2, 1, 1) \cdot (-1, 2, 0) = -2 + 2 = 0$. Коэффициенты $(y_1 - y_2,\ x_2 - x_1,\ x_1 y_2 - x_2 y_1)$ прямой через две точки те же, что в разделе 4.1.3 «Уравнения прямой».

Точка пересечения. Пусть $l$ и $m$ — прямые и $S = l \times m$. По тождествам $S \cdot l = 0$ и $S \cdot m = 0$ вектор $S$, если он ненулевой, лежит на обеих прямых. Если $l$ и $m$ не пропорциональны, то $S \neq 0$ по факту, и общие точки прямых — ненулевые $Z$ с $Z \cdot l = 0$ и $Z \cdot m = 0$. Два независимых линейных уравнения от трёх неизвестных задают одномерное подпространство, натянутое на $S$; значит, общая точка одна с точностью до множителя: $S$.

Координата $W$ точки $S$ равна $l_1 m_2 - l_2 m_1$. Если $W = 0$, то $(l_1, l_2)$ и $(m_1, m_2)$ пропорциональны (по тому же рассуждению, что в факте, при $(l_1, l_2) \neq (0, 0)$), и прямые параллельны или совпадают. Параллельные прямые ($S \neq 0$) имеют общую точку на бесконечности $S = (S_1, S_2, 0)$. Если $S = 0$, то $m = \mu l$ по факту, и прямые совпадают: тогда общих точек бесконечно много, и $S$ их не задаёт.

Пример. Пусть $l = P \times Q = (-1, 2, 0)$ из предыдущего примера и $m = (1, 0, -1)$ — прямая $x = 1$. Тогда $$S = l \times m = (2 \cdot (-1) - 0 \cdot 0,\; 0 \cdot 1 - (-1)(-1),\; (-1) \cdot 0 - 2 \cdot 1) = (-2, -1, -2).$$ Каноническая тройка получается умножением на $-1$: $(2, 1, 2)$, здесь $W = 2 > 0$ и $\operatorname{gcd}(2, 1, 2) = 1$. Значит, точка пересечения имеет координаты $x = 2/2 = 1$ и $y = 1/2$, то есть $(1, 1/2)$. Проверка: $(-1, 2, 0) \cdot (2, 1, 2) = -2 + 2 = 0$ и $(1, 0, -1) \cdot (2, 1, 2) = 2 - 2 = 0$. Тот же результат дают формулы Крамера из раздела «Пересечение 2-х прямых»; векторное произведение здесь удобнее, потому что случай параллельных прямых виден сразу: $W = 0$.

Двойственность и ориентация

Точки и прямые — объекты одного рода: и то и другое задаётся тройкой, а отношение «точка лежит на прямой» записывается одной формулой $P \cdot l = 0$. Поэтому одна и та же операция $\times$ строит прямую по двум точкам и точку по двум прямым.

Для точек $P = (x_1, y_1, 1)$, $Q = (x_2, y_2, 1)$, $R = (x_3, y_3, 1)$ обозначим через $\det(P, Q, R)$ определитель матрицы, строками которой служат $P$, $Q$, $R$. Три точки лежат на одной прямой тогда и только тогда, когда $\det(P, Q, R) = 0$.

Доказательство. Разложив определитель по строке $R$, получаем $\det(P, Q, R) = (P \times Q) \cdot R$. Если $P = Q$, то две строки совпадают, определитель равен нулю, и точки очевидно лежат на одной прямой. Если $P \neq Q$, то прямая $P \times Q$ — единственная с точностью до множителя прямая через $P$ и $Q$ (упражнение 3), и $R$ лежит на ней тогда и только тогда, когда $(P \times Q) \cdot R = 0$.

Для плоских векторов $s = (s_1, s_2)$ и $t = (t_1, t_2)$ положим $s \times t = s_1 t_2 - s_2 t_1$; это третья компонента векторного произведения $(s_1, s_2, 0) \times (t_1, t_2, 0)$. Вычтем первую строку из второй и третьей: определитель не изменится. В последнем столбце останутся $1, 0, 0$, и разложение по нему даёт $$\det(P, Q, R) = (x_2 - x_1)(y_3 - y_1) - (y_2 - y_1)(x_3 - x_1) = (Q - P) \times (R - P).$$

Пример. Для $A = (1, 1)$, $B = (3, 2)$, $C = (5, 3)$ имеем $(B - A) \times (C - A) = 2 \cdot 2 - 1 \cdot 4 = 0$; действительно, $C - A = 2(B - A)$, и точки лежат на одной прямой. Для $A = (0, 0)$, $B = (4, 0)$, $C = (3, 3)$ имеем $(B - A) \times (C - A) = 4 \cdot 3 - 0 \cdot 3 = 12$, и то же число получается из прямой через $A$ и $B$: $(A \times B) \cdot C = (0, 4, 0) \cdot (3, 3, 1) = 12$.

Знак величины $(Q - P) \times (R - P)$ задаёт поворот. Поскольку $R - Q = (R - P) - (Q - P)$, а $(Q - P) \times (Q - P) = 0$, то $$(Q - P) \times (R - Q) = (Q - P) \times (R - P),$$ и знак показывает поворот пути $P \to Q \to R$ в точке $Q$: больше нуля — поворот влево (против часовой стрелки), меньше нуля — вправо, ноль — точки лежат на одной прямой. Почему знак отвечает за поворот: если $s = |s|(\cos\alpha, \sin\alpha)$ и $t = |t|(\cos\beta, \sin\beta)$, то $$\begin{aligned} s \times t &= s_1 t_2 - s_2 t_1 \\ &= |s|\,|t|\,(\cos\alpha \sin\beta - \sin\alpha \cos\beta) \\ &= |s|\,|t|\sin(\beta - \alpha), \end{aligned}$$ и знак совпадает со знаком $\sin(\beta - \alpha)$. Если $0 < \beta - \alpha < \pi$, то вектор $t$ получается из $s$ поворотом против часовой стрелки на угол меньше развёрнутого, то есть поворотом влево. Например, для $s = (1, 0)$ и $t = (0, 1)$ имеем $s \times t = 1 > 0$: от направления «на восток» к направлению «на север» поворачиваем влево. В лекции об алгоритмах Грэхема и Джарвиса, где строится выпуклая оболочка, эта же величина записана как $(a - o) \times (b - o)$ при $o = P$, $a = Q$, $b = R$, и служит проверкой левого поворота. Ориентация тройки точек разобрана также в разделе «Ориентация тройки точек».

Вычисления без ошибок округления

Если координаты входных точек и прямых целые, то все величины этой лекции — векторные и скалярные произведения — целые, поскольку получаются только сложением и умножением. Нормировка точки пересечения использует лишь точное деление на общий делитель, поэтому дробей внутри программ не возникает; дробная запись появляется только при печати.

Тройка $(X, Y, W)$ с целыми координатами и $W \neq 0$ называется канонической, если $W > 0$ и $\operatorname{gcd}(X, Y, W) = 1$, где $\operatorname{gcd}(X, Y, W)$ — наибольший общий делитель модулей $X$, $Y$, $W$.

Утверждение. Точка с $W \neq 0$, заданная целой тройкой, имеет ровно одну каноническую тройку.

Доказательство. Существование: делим тройку на $g = \operatorname{gcd}(X, Y, W)$, а при $W < 0$ берём $-g$; получаем тройку с положительным $W$ и $\operatorname{gcd} = 1$. Единственность: пусть $v$ и $v'$ — канонические тройки одной точки. Тогда $v' = \lambda v$ с $\lambda = W'/W > 0$. Запишем $\lambda = p/q$ с натуральными взаимно простыми $p$ и $q$. Из $q v' = p v$ следует, что $q$ делит все координаты $v$ (так как $q$ взаимно просто с $p$), а значит, делит $\operatorname{gcd}(v) = 1$, то есть $q = 1$. Теперь $v' = p v$, и $p$ делит все координаты $v'$, поэтому $p = 1$. Итак, $v = v'$.

Следствие. Две точки с $W \neq 0$ совпадают тогда и только тогда, когда совпадают их канонические тройки. Функция intersection возвращает каноническую тройку, поэтому равенство точек пересечения проверяется сравнением троек.

Кроме нормировки, каждая функция выполняет фиксированное число арифметических операций: векторное произведение — шесть умножений и три вычитания. Нормировка выполняет два вычисления $\operatorname{gcd}$ и три деления; алгоритм Евклида для $\operatorname{gcd}$ делает $O(\log M)$ делений, где $M$ — наибольший по модулю вход. При фиксированном $M$ это константа.

Граница 64-битной арифметики. Программы на C++ и Rust используют 64-битные знаковые целые: диапазон от $-2^{63}$ до $2^{63} - 1$, где $2^{63} \approx 9{,}22 \cdot 10^{18}$. Программа на Python использует целые числа произвольной длины и точна при любых входах.

Пусть каждая входная точка имеет вид $(x, y, 1)$ с $|x| \le 10^6$ и $|y| \le 10^6$. Тогда прямая $l = P \times Q$ имеет $|l_1|, |l_2| \le 2 \cdot 10^6$ и $|l_3| \le 2 \cdot 10^{12}$. Компоненты $S_1 = l_2 m_3 - l_3 m_2$ и $S_2 = l_3 m_1 - l_1 m_3$ — разности двух произведений, в каждом из которых один множитель имеет индекс 3 (не больше $2 \cdot 10^{12}$), а другой — индекс 1 или 2 (не больше $2 \cdot 10^6$); каждое такое произведение по модулю не больше $4 \cdot 10^{18}$. Компонента $S_3 = l_1 m_2 - l_2 m_1$ — разность произведений индексов 1 и 2, каждое не больше $4 \cdot 10^{12}$. Поэтому $|S_1|, |S_2| \le 8 \cdot 10^{18}$, а $|S_3| \le 8 \cdot 10^{12}$. Так как $8 \cdot 10^{18} < 2^{63}$, все промежуточные значения line_through и intersection помещаются в 64 бита, и цепочка «прямая через две точки, затем пересечение» точна. Для входной точки $P$ значения $P \cdot l$, где $l$ — прямая через входные точки, не превосходят $6 \cdot 10^{12}$; для цепочки «точка пересечения, затем on_line» эта оценка не действует (см. ниже).

При $|x|, |y| \le 10^6$ координаты точки пересечения не больше $8 \cdot 10^{18} < 2^{63}$, и они помещаются в 64 бита. Но 64 бит не хватает в двух случаях: для значений, которые подставляются в дальнейшие вычисления (цепочка «точка пересечения, затем on_line», см. ниже), и при $|x|, |y| \ge 2^{20}$, когда сама точка пересечения может иметь координату не меньше $2^{63}$. Тогда нужны 128-битные целые (__int128 в C++, i128 в Rust) или Python. Например, цепочка «точка пересечения, затем on_line» при $10^6$ уже может выйти за пределы 64 бит: её промежуточные значения не превосходят $48 \cdot 10^{24}$. Для этой цепочки 64-битная арифметика точна при $|x|, |y| \le 2 \cdot 10^4$: при $B = 2 \cdot 10^4$ координаты точки пересечения не превосходят $8B^3$, коэффициенты прямой — $2B$ и $2B^2$, и все промежуточные значения не больше $48B^4 = 7{,}68 \cdot 10^{18} < 2^{63}$.

Упражнения

1. Найдите прямую, проходящую через точку $(1, 2)$ параллельно прямой $l = (3, -1, 5)$.

Ответ: $(3, -1, -1)$, то есть $3x - y - 1 = 0$. Проверка: у ответа те же $(a, b) = (3, -1)$, что у $l$, поэтому прямые параллельны; свободные члены $-1$ и $5$ различны, значит, прямые различны. Далее $(1, 2, 1) \cdot (3, -1, -1) = 3 - 2 - 1 = 0$, то есть точка $(1, 2)$ лежит на ответе. Другой способ: точка на бесконечности прямой $l$ равна $D = (1, 3, 0)$, так как $(1, 3, 0) \cdot l = 3 - 3 = 0$; тогда $$(1, 2, 1) \times (1, 3, 0) = (-3, 1, 1) = -(3, -1, -1).$$

2. Найдите точку пересечения диагоналей четырёхугольника с вершинами $A = (0, 0)$, $B = (4, 0)$, $C = (3, 3)$, $D = (0, 2)$, перечисленными в порядке обхода.

Ответ: $(4/3, 4/3)$. Диагонали — $AC$ и $BD$. Прямая $AC$: $$(0, 0, 1) \times (3, 3, 1) = (-3, 3, 0) = -3(1, -1, 0),$$ то есть $x - y = 0$. Прямая $BD$: $$(4, 0, 1) \times (0, 2, 1) = (-2, -4, 8) = -2(1, 2, -4),$$ то есть $x + 2y - 4 = 0$. Их пересечение: $(1, -1, 0) \times (1, 2, -4) = (4, 4, 3)$, то есть точка $(4/3, 4/3)$. Проверка: $4/3 - 4/3 = 0$ и $4/3 + 8/3 - 4 = 0$.

3. Докажите, что при $P \neq Q$ прямая через $P$ и $Q$ определена однозначно с точностью до ненулевого множителя.

Ответ. Пусть $l$ и $l'$ — ненулевые векторы прямых через $P$ и $Q$, то есть $l \cdot P = l \cdot Q = 0$ и $l' \cdot P = l' \cdot Q = 0$. Векторы $P$ и $Q$ линейно независимы: иначе $Q = \mu P$ при некотором $\mu$, а так как $Q \neq 0$, то $\mu \neq 0$, то есть $P = Q$ как точки. Множество $\{Z : Z \cdot P = 0,\ Z \cdot Q = 0\}$ — подпространство пространства троек, заданное двумя независимыми линейными уравнениями, и его размерность равна $3 - 2 = 1$. Векторы $l$ и $l'$ ненулевые и лежат в этом одномерном подпространстве, поэтому $l' = \mu l$ при некотором $\mu \neq 0$. $\blacksquare$

Программы

Ниже приведены три полные программы: на Python, C++ и Rust. Они реализуют одни и те же функции и печатают одинаковую демонстрацию: прямую через точки $(0, 0)$ и $(2, 1)$, точку пересечения с прямой $x = 1$ и результат проверки on_line для обеих точек. Функции:

• cross3(u, v) — векторное произведение $u \times v$;
• dot3(u, v) — скалярное произведение $u \cdot v$;
• point(x, y) — точка $(x, y, 1)$;
• line_through(p, q) — прямая $p \times q$ (предполагается $p \neq q$ и хотя бы одна точка конечна);
• intersection(l, m) — каноническая тройка точки пересечения; если прямые параллельны или совпадают, функция не возвращает значения: это None в Python и Rust (в Rust тип Option) и std::nullopt в C++;
• on_line(p, l) — проверка $p \cdot l = 0$;
• point_text(p) — текстовый вид точки с $W \neq 0$, например (1, 1/2).

Вспомогательные функции: normalize (каноническая тройка), fraction_text (дробь), triple_text (тройка целых чисел). В C++ и Rust есть ещё bool_text (текст логического значения), а в Rust — gcd (наибольший общий делитель). Все входные числа демонстрации по модулю не больше 2, то есть далеко внутри границы из раздела «Вычисления без ошибок округления».

Реализация на Python

# Однородные координаты на плоскости: точки, прямые, пересечения.
# Только целые числа (произвольной длины) и стандартная библиотека. Запуск: python3 -I main.py
from math import gcd


def cross3(u, v):
    """Векторное произведение u x v троек (u1, u2, u3) и (v1, v2, v3)."""
    u1, u2, u3 = u
    v1, v2, v3 = v
    return (u2 * v3 - u3 * v2, u3 * v1 - u1 * v3, u1 * v2 - u2 * v1)


def dot3(u, v):
    """Скалярное произведение u . v."""
    return u[0] * v[0] + u[1] * v[1] + u[2] * v[2]


def point(x, y):
    """Точка плоскости (x, y) в однородных координатах: (x, y, 1)."""
    return (x, y, 1)


def line_through(p, q):
    """Прямая через точки p и q (при p != q, хотя бы одна точка конечна): вектор p x q."""
    return cross3(p, q)


def normalize(p):
    """Каноническая тройка точки с W != 0: W > 0 и gcd(X, Y, W) = 1."""
    x, y, w = p
    g = gcd(gcd(x, y), w)
    if w < 0:
        g = -g
    return (x // g, y // g, w // g)


def intersection(l, m):
    """Точка пересечения прямых l и m или None, если они параллельны или совпадают."""
    p = cross3(l, m)
    if p[2] == 0:
        return None
    return normalize(p)


def on_line(p, l):
    """Точка p лежит на прямой l тогда и только тогда, когда p . l = 0."""
    return dot3(p, l) == 0


def fraction_text(num, den):
    """Дробь num/den (den != 0) в виде '3' или '-1/2'."""
    if den < 0:
        num, den = -num, -den
    g = gcd(num, den)
    num, den = num // g, den // g
    return str(num) if den == 1 else f"{num}/{den}"


def point_text(p):
    """Текстовый вид точки с W != 0, например (1, 1/2)."""
    x, y, w = p
    return f"({fraction_text(x, w)}, {fraction_text(y, w)})"


def triple_text(v):
    """Тройка целых чисел в виде (a, b, c)."""
    return "({}, {}, {})".format(*v)


def main():
    p = point(0, 0)
    q = point(2, 1)
    l = line_through(p, q)
    m = (1, 0, -1)  # прямая x = 1
    x = intersection(l, m)
    print(f"line through (0, 0) and (2, 1): {triple_text(l)}")
    print("intersection with x = 1: " + ("none" if x is None else point_text(x)))
    print(f"on_line((0, 0), l): {on_line(p, l)}")
    print(f"on_line((2, 1), l): {on_line(q, l)}")


if __name__ == "__main__":
    main()

Реализация на C++

// Однородные координаты на плоскости: точки, прямые, пересечения.
// Только целые числа (int64): для точек (x, y, 1) с |x|, |y| <= 10^6 все значения помещаются
// в 64 бита, для цепочек см. текст лекции. Сборка: g++ -std=c++17 -Wall -Wextra -Werror -O2 main.cpp
#include <cstdint>
#include <iostream>
#include <numeric>
#include <optional>
#include <string>

using i64 = std::int64_t;

// Тройка (x, y, z): точка (X, Y, W) или прямая (a, b, c); z — это W или c
struct Vec3 {
    i64 x;
    i64 y;
    i64 z;
};

// Векторное произведение u x v
Vec3 cross3(Vec3 u, Vec3 v) {
    return {u.y * v.z - u.z * v.y, u.z * v.x - u.x * v.z, u.x * v.y - u.y * v.x};
}

// Скалярное произведение u . v
i64 dot3(Vec3 u, Vec3 v) { return u.x * v.x + u.y * v.y + u.z * v.z; }

// Точка плоскости (x, y) в однородных координатах: (x, y, 1)
Vec3 point(i64 x, i64 y) { return {x, y, 1}; }

// Прямая через точки p и q (при p != q, хотя бы одна точка конечна): вектор p x q
Vec3 line_through(Vec3 p, Vec3 q) { return cross3(p, q); }

// Каноническая тройка точки с W != 0: W > 0 и gcd(X, Y, W) = 1
Vec3 normalize(Vec3 p) {
    const i64 g0 = std::gcd(std::gcd(p.x, p.y), p.z);
    const i64 g = p.z < 0 ? -g0 : g0;
    return {p.x / g, p.y / g, p.z / g};
}

// Точка пересечения прямых l и m или std::nullopt, если они параллельны или совпадают
std::optional<Vec3> intersection(Vec3 l, Vec3 m) {
    const Vec3 p = cross3(l, m);
    if (p.z == 0) return std::nullopt;
    return normalize(p);
}

// Точка p лежит на прямой l тогда и только тогда, когда p . l = 0
bool on_line(Vec3 p, Vec3 l) { return dot3(p, l) == 0; }

// Дробь num/den (den != 0) в виде "3" или "-1/2"
std::string fraction_text(i64 num, i64 den) {
    if (den < 0) {
        num = -num;
        den = -den;
    }
    const i64 g = std::gcd(num, den);
    num /= g;
    den /= g;
    return den == 1 ? std::to_string(num) : std::to_string(num) + "/" + std::to_string(den);
}

// Текстовый вид точки с W != 0, например (1, 1/2)
std::string point_text(Vec3 p) {
    return "(" + fraction_text(p.x, p.z) + ", " + fraction_text(p.y, p.z) + ")";
}

// Тройка целых чисел в виде (a, b, c)
std::string triple_text(Vec3 v) {
    return "(" + std::to_string(v.x) + ", " + std::to_string(v.y) + ", " + std::to_string(v.z) +
           ")";
}

const char* bool_text(bool b) { return b ? "True" : "False"; }

int main() {
    const Vec3 p = point(0, 0);
    const Vec3 q = point(2, 1);
    const Vec3 l = line_through(p, q);
    const Vec3 m = {1, 0, -1};  // прямая x = 1
    const std::optional<Vec3> x = intersection(l, m);
    std::cout << "line through (0, 0) and (2, 1): " << triple_text(l) << "\n";
    std::cout << "intersection with x = 1: " << (x ? point_text(*x) : std::string("none")) << "\n";
    std::cout << "on_line((0, 0), l): " << bool_text(on_line(p, l)) << "\n";
    std::cout << "on_line((2, 1), l): " << bool_text(on_line(q, l)) << "\n";
    return 0;
}

Реализация на Rust

// Однородные координаты на плоскости: точки, прямые, пересечения.
// Только целые числа (i64): для точек (x, y, 1) с |x|, |y| <= 10^6 все значения помещаются
// в 64 бита, для цепочек см. текст лекции. Сборка: rustc --edition 2021 -D warnings -O main.rs

// Тройка (x, y, z): точка (X, Y, W) или прямая (a, b, c); z — это W или c
pub type Vec3 = (i64, i64, i64);

// Векторное произведение u x v
pub fn cross3(u: Vec3, v: Vec3) -> Vec3 {
    let (u1, u2, u3) = u;
    let (v1, v2, v3) = v;
    (u2 * v3 - u3 * v2, u3 * v1 - u1 * v3, u1 * v2 - u2 * v1)
}

// Скалярное произведение u . v
pub fn dot3(u: Vec3, v: Vec3) -> i64 {
    u.0 * v.0 + u.1 * v.1 + u.2 * v.2
}

// Точка плоскости (x, y) в однородных координатах: (x, y, 1)
pub fn point(x: i64, y: i64) -> Vec3 {
    (x, y, 1)
}

// Прямая через точки p и q (при p != q, хотя бы одна точка конечна): вектор p x q
pub fn line_through(p: Vec3, q: Vec3) -> Vec3 {
    cross3(p, q)
}

fn gcd(a: i64, b: i64) -> i64 {
    let (mut a, mut b) = (a.abs(), b.abs());
    while b != 0 {
        (a, b) = (b, a % b);
    }
    a
}

// Каноническая тройка точки с W != 0: W > 0 и gcd(X, Y, W) = 1
fn normalize(p: Vec3) -> Vec3 {
    let g = gcd(gcd(p.0, p.1), p.2);
    let g = if p.2 < 0 { -g } else { g };
    (p.0 / g, p.1 / g, p.2 / g)
}

// Точка пересечения прямых l и m или None, если они параллельны или совпадают
pub fn intersection(l: Vec3, m: Vec3) -> Option<Vec3> {
    let p = cross3(l, m);
    if p.2 == 0 {
        return None;
    }
    Some(normalize(p))
}

// Точка p лежит на прямой l тогда и только тогда, когда p . l = 0
pub fn on_line(p: Vec3, l: Vec3) -> bool {
    dot3(p, l) == 0
}

// Дробь num/den (den != 0) в виде "3" или "-1/2"
fn fraction_text(num: i64, den: i64) -> String {
    let (num, den) = if den < 0 { (-num, -den) } else { (num, den) };
    let g = gcd(num, den);
    let (num, den) = (num / g, den / g);
    if den == 1 {
        num.to_string()
    } else {
        format!("{num}/{den}")
    }
}

// Текстовый вид точки с W != 0, например (1, 1/2)
pub fn point_text(p: Vec3) -> String {
    format!("({}, {})", fraction_text(p.0, p.2), fraction_text(p.1, p.2))
}

// Тройка целых чисел в виде (a, b, c)
pub fn triple_text(v: Vec3) -> String {
    format!("({}, {}, {})", v.0, v.1, v.2)
}

pub fn bool_text(b: bool) -> &'static str {
    if b {
        "True"
    } else {
        "False"
    }
}

fn main() {
    let p = point(0, 0);
    let q = point(2, 1);
    let l = line_through(p, q);
    let m = (1, 0, -1); // прямая x = 1
    println!("line through (0, 0) and (2, 1): {}", triple_text(l));
    match intersection(l, m) {
        Some(x) => println!("intersection with x = 1: {}", point_text(x)),
        None => println!("intersection with x = 1: none"),
    }
    println!("on_line((0, 0), l): {}", bool_text(on_line(p, l)));
    println!("on_line((2, 1), l): {}", bool_text(on_line(q, l)));
}

Все три программы печатают один и тот же результат. Его числа совпадают с разделом «Прямая через две точки и точка пересечения»: прямая $(-1, 2, 0)$, точка $(1, 1/2)$, и обе точки $(0, 0)$ и $(2, 1)$ лежат на прямой.

line through (0, 0) and (2, 1): (-1, 2, 0)
intersection with x = 1: (1, 1/2)
on_line((0, 0), l): True
on_line((2, 1), l): True

Преобразования плоскости — сдвиг, поворот, растяжение и проективные преобразования — записываются матрицами $3 \times 3$ и описаны в следующей лекции «Преобразования плоскости в однородных координатах».