Точность: когда float врёт о повороте и что с этим делать

Если ты знаешь, почему поворот в float ошибается около прямой, в каком диапазоне целые координаты дают точный ответ и почему eps-сравнение не лечит, листай до «Теперь сам».

Как думаешь, сколько ошибок

Возьмём прямую y = x и две точки на ней, (12, 12) и (24, 24). Около точки (0,5; 0,5) выложим сетку 16 × 16 точек с самым маленьким шагом, какой бывает у float в этом месте, - 2⁻⁵³ ≈ 1,1 · 10⁻¹⁶ (у float 53 значащих бита, и у чисел от 0,5 до 1 младший весит ровно 2⁻⁵³). Для каждой точки спросим знак поворота относительно прямой. Правильный ответ очевиден: точки выше диагонали - слева, ниже - справа, на диагонали - ноль. Картина - два треугольника и диагональ нулей.

Ответь себе до чтения: сколько из 256 знаков float посчитает неверно? Ни одного, потому что двойная точность - это 16 знаков? Пару штук на самой диагонали? Ответ: 114 - почти половину. И ошибки не случайный шум у диагонали, а целые полосы неверных ответов.

На что это похоже

Плотник проверяет, стоят ли три гвоздя на одной линии, линейкой с миллиметровыми делениями. Если гвоздь отклонился на два сантиметра, линейка это покажет. Если на десятую миллиметра - ответ зависит от того, как легла линейка, и два плотника скажут разное. Хуже того, один плотник, проверив гвозди парами в разном порядке, сам себе противоречит: «первый, второй и третий на линии», «второй, третий и четвёртый на линии», а первый, второй и четвёртый - нет.

float - это такая линейка с очень мелкими делениями, но всё-таки с делениями. А геометрический алгоритм - это плотник, который строит дом по своим замерам. Противоречивые замеры не дают «немного кривой» дом: они дают дом, у которого стена проходит сквозь крышу.

Механизм: предикаты, конструкции и вырожденность

Операции геометрии делятся на два вида. Предикат отвечает на вопрос с дискретным ответом: знак поворота, пересекаются ли отрезки, внутри ли точка. От него зависит ветвление программы. Конструкция строит новый объект: точку пересечения, центр окружности, расстояние. Ошибка в конструкции - это точка, нарисованная на крупинку не там. Ошибка в предикате - это другая ветка if, и алгоритм дальше работает с неверной картиной мира.

Почему float ошибается именно в предикатах. В поворот входят произведения: здесь они порядка 11,5 · 23,5 ≈ 270. У чисел такого размера шаг float около 5,7 · 10⁻¹⁴, и каждое произведение округляется до этого шага. А настоящая величина поворота здесь 12 · (b − a) - для соседних точек сетки порядка 10⁻¹⁵, в десятки раз меньше погрешности. Результат - разность двух округлённых чисел, и знак этой разности решает округление, а не геометрия.

Такие случаи - три точки на одной прямой, четыре на одной окружности, два одинаковых x - называют вырожденными. Учебники часто описывают алгоритм в предположении общего положения: вырожденных для этого алгоритма случаев нет - чаще всего требуют, чтобы никакие три точки не лежали на одной прямой. Реальные данные этому не подчиняются - сетки, здания, выровненные по улице, точки, введённые вручную с шагом в метр. Значит, алгоритм обязан и правильно обрабатывать вырожденные случаи, и узнавать их, а узнать их можно только точным предикатом.

ВХОД                        ПОВОРОТ НЕ БОЛЬШЕ          ТОЧЕН?
целые, |x|, |y| <= 2^25     8 · 2^50 = 2^53            да, даже в float
целые, |x|, |y| < 2^30      меньше 2^63                да, в int64
целые любой величины        сколько угодно             да, в целых Python
десятичные дроби            -                          нет, если в float

Граница в первой строке считается так: координаты по модулю не больше M, разности не больше 2M, каждое произведение не больше 4M², разность двух произведений не больше 8M². Целые числа до 2⁵³ float хранит точно, поэтому при 8M² ≤ 2⁵³, то есть M ≤ 2²⁵ ≈ 3,3 · 10⁷, поворот целых координат не ошибается даже в float.

Правило

Решения - знак поворота, пересечение, «внутри» - принимай только точными предикатами: целые координаты в безопасном диапазоне (до 2²⁵ для float, меньше 2³⁰ для 64-битных целых), рациональные числа или адаптивные предикаты. Приближённо можно считать только конструкции, по которым ничего не решают. Допуск eps - свойство данных: округли координаты к сетке один раз на входе (шаг сетки - точность самих измерений), а дальше считай точно.

Разбор: картина знаков около прямой

Сетка из «Как думаешь». Точный знак здесь считается без огромных чисел: для p = (a, b) поворот (q − p) × (r − p) при q = (12, 12), r = (24, 24) раскрывается в 12 · (b − a), то есть знак j − i.

def povorot(a, b, c):
    return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])

def simvol(z):
    if z > 0:
        return "+"
    if z < 0:
        return "-"
    return "0"

u = 2 ** -53
q = [12.0, 12.0]
r = [24.0, 24.0]
oshibok = 0
for j in range(15, -1, -1):
    stroka = ""
    tochno = ""
    for i in range(16):
        p = [0.5 + i * u, 0.5 + j * u]
        s = simvol(povorot(p, q, r))
        t = simvol(j - i)
        stroka = stroka + s
        tochno = tochno + t
        if s != t:
            oshibok = oshibok + 1
    print(stroka, tochno)
print(oshibok)

+++++++++0000000 +++++++++++++++0
+++++++++0000000 ++++++++++++++0-
+++++++++0000000 +++++++++++++0--
+++++++++0000000 ++++++++++++0---
+++++++++0000000 +++++++++++0----
+++++++++0000000 ++++++++++0-----
+++++++++0000000 +++++++++0------
000000000------- ++++++++0-------
000000000------- +++++++0--------
000000000------- ++++++0---------
000000000------- +++++0----------
000000000------- ++++0-----------
000000000------- +++0------------
000000000------- ++0-------------
000000000------- +0--------------
000000000------- 0---------------
114

Слева - что посчитал float, справа - правда, внизу - число расхождений. Вместо диагонали - четыре прямоугольных блока. Левый верхний и правый нижний верны: там точка далеко от прямой по меркам сетки. А два других блока целиком из нулей: квадраты 7 × 7 и 9 × 9 точек «лежат на прямой», хотя на прямой из них только диагональ. Это 42 и 72 ошибки - те самые 114. Хуже всего не число ошибок, а их вид: множество точек «на прямой» оказалось квадратом, а прямая квадратом быть не может. Алгоритм оболочки, получивший такие ответы, может выкинуть точку, лежащую снаружи, или зациклиться, переставляя две точки туда-обратно.

Важно и то, чего здесь нет: никаких экзотических чисел. Это самые обычные float, сетка лишь мелкая. В реальных данных та же ситуация - три точки, почти лежащие на прямой, - возникает, когда к точкам пересечения, уже округлённым, добавляют новые проверки.

Разбор: допуск eps не транзитивен

Первая мысль - «раз точно нельзя, сравню с нулём не точно»: считать, что три точки на прямой, если точка отстоит от хорды не больше чем на eps. Ломаная ниже - плавная дуга: каждая точка отстоит от хорды через своих соседей меньше чем на 0,5, а от хорды через концы - почти на 3. «Почти на прямой» для соседей не значит «почти на прямой» для всех, то есть отношение не транзитивно. Посмотрим, что это делает с обычной операцией - упрощением ломаной, которое выкидывает «лишние» средние точки.

import math

def povorot(a, b, c):
    return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])

def pochti_pryamo(a, b, c, eps):
    # расстояние от b до прямой ac не больше eps
    dlina = math.sqrt((c[0] - a[0]) ** 2 + (c[1] - a[1]) ** 2)
    return abs(povorot(a, b, c)) <= eps * dlina

def uprostit(tochki, eps):
    itog = [tochki[0]]
    for p in tochki[1:]:
        if len(itog) >= 2 and pochti_pryamo(itog[-2], itog[-1], p, eps):
            itog.pop()
        itog.append(p)
    return itog

lomanaya = [[0, 0], [10, 0], [20, 1], [30, 3], [40, 6], [50, 10]]
for eps in [0, 1]:
    vpered = uprostit(lomanaya, eps)
    nazad = uprostit(lomanaya[::-1], eps)[::-1]
    print(eps, vpered)
    print(eps, nazad)

0 [[0, 0], [10, 0], [20, 1], [30, 3], [40, 6], [50, 10]]
0 [[0, 0], [10, 0], [20, 1], [30, 3], [40, 6], [50, 10]]
1 [[0, 0], [30, 3], [50, 10]]
1 [[0, 0], [20, 1], [50, 10]]

С точным сравнением (eps = 0) обе стороны обхода дают одно и то же: ни одна точка не лежит на прямой своих соседей, выкидывать нечего. С eps = 1 проход слева направо оставил излом в (30, 3), справа налево - в (20, 1). Оба ответа законны: каждая выкинутая точка не дальше 1 от своей хорды. Но они разные.

Где это бьёт. Два соседних участка на карте делят общую границу. В контуре одного участка граница записана в одну сторону, в контуре другого - в обратную: так устроен обход против часовой стрелки у обоих. Упрости каждый участок отдельно - и общая граница превратится в две разные линии, а между участками появятся щели и наложения. Поэтому в библиотеках есть отдельное упрощение покрытий, которое упрощает общую границу один раз для обоих соседей: ST_CoverageSimplify в PostGIS 3.4, coverage_simplify в Shapely 2.1.

Что ломается без этого

abs(z) < 1e-9 вместо z == 0. Порог не масштабируется: для координат в километрах он слишком мелкий, для координат в градусах - огромный. И он ломает транзитивность, как во втором разборе.

Округление после вычислений. Точки пересечения округлили до сантиметра, чтобы сохранить в базу, - и многоугольник стал самопересекающимся: вершина переехала через соседнее ребро. Правильный порядок обратный: сначала сетка, потом точная геометрия на ней. Есть и специальная техника - округление к сетке с сохранением топологии (snap rounding): она двигает вершины к узлам сетки так, что рёбра не перескакивают друг через друга.

Переполнение там, где «целые - значит точные». Координаты OpenStreetMap хранятся целыми в десятимиллионных долях градуса: долгота от −1,8 · 10⁹ до 1,8 · 10⁹. Это больше 2³⁰, и поворот по таким целым на всём земном шаре в 64 бита уже не помещается. В пределах одного города разности маленькие, и поворот, который умножает только разности, точен. А формула шнурка умножает сами координаты - для неё сначала вычитают опорную точку, а потом умножают.

Зачем это в работе

Ты выбираешь, как хранить координаты, - это выбор точности всех будущих проверок. Целые в фиксированных единицах (миллиметры, 10⁻⁷ градуса) дают точные предикаты почти даром. Если координаты дробные, в Python их переводят в целые округлением round(x * 10 ** k) - лучше прямо из десятичной строки, пока точность не потеряна, - или считают в fractions.Fraction, а большую геометрию отдают библиотеке, которая знает об этой проблеме: в CGAL есть ядро Exact_predicates_inexact_constructions_kernel - ровно правило этого урока, записанное в имени типа. Предикаты Шевчука (1996) делают то же адаптивно: сначала считают в float и оценивают погрешность, а точную арифметику включают, только если результат меньше погрешности. На типичных данных это редкость; на данных с сеткой и точками на одной прямой - чаще, и тогда расчёт дороже, но ответ остаётся верным.

Теперь сам

1. Координаты складов - целые от −10⁶ до 10⁶. Точен ли поворот в float?
2. Что из этого предикат, а что конструкция: точка пересечения двух отрезков; «пересекаются ли»; площадь многоугольника; «по часовой ли обход»?
3. Почему формулу шнурка на координатах OpenStreetMap безопаснее считать от первой вершины, а не от начала координат?