как вы можете практически решить выпуклую оболочку, когда есть проблемы с плавающей точностью?

Предположим, что у вас есть 100000 точек на кривой y = x^2. Вы хотите найти выпуклую оболочку этих точек. Все координаты являются плавающими числами.

В моей реализации сканирования Грэма единственное место, где я работаю с плавающими числами, — это когда я сначала сортирую все точки по их координатам, а затем у меня есть одна функция, которая определяет, делают ли три точки поворот влево или вправо.

Точки:

struct point {
   double x; 
   double y;
};

Сортировочный компаратор:

inline bool operator() (const point &p1, const point &p2) {
    return (p1.x < p2.x) || (p1.x == p2.x && p1.y > p2.y);
}

Левый/правый поворот:

inline int ccw(point *p1, point *p2, point *p3) {
  double left  = (p1->x - p3->x)*(p2->y - p3->y);
  double right = (p1->y - p3->y)*(p2->x - p3->x);
  double res = left - right;
  return res > 0;
}

Моя программа говорит, что из 100 000 точек только 68894 являются частью выпуклой оболочки. Но поскольку они находятся на кривой, все они должны быть частью выпуклой оболочки.

Для вашего глаза это не имеет никакого значения. См. рисунок ниже. Красные точки являются частью выпуклой оболочки.

изображение

Но если вы посмотрите достаточно близко и приблизите точки, вы увидите, что некоторые из них синие, поэтому они не включены в выпуклую оболочку.

изображение

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

Думаю, я мог бы использовать внешнюю библиотеку с произвольной точностью для чисел с плавающей запятой, но меня больше интересуют простые типы данных, которые есть, например, в C++.

Как я могу повысить точность? Я читал об эпсилонах, но как здесь поможет использование эпсилонов? Я бы по-прежнему считал некоторые точки, которые близки друг к другу, одинаковыми, поэтому я не получу точность ближе к 100%.

Как лучше всего подойти к этой проблеме?


person jsguy    schedule 30.09.2014    source источник
comment
ты пробовал long double?   -  person mch    schedule 30.09.2014
comment
Возможно ли, что ваш алгоритм выпуклой оболочки верен, но это округление произошло, когда вы первоначально оценивали y = x ^ 2?   -  person ajclinto    schedule 30.09.2014
comment
Во-первых, никакая конечная точность не позволит вам представлять действительные числа. Во-вторых, точки, которые вы рисуете, являются приблизительными. В-третьих, настоящие конструктивные реалы (бесконечная точность) нельзя сравнивать так, как вы хотите их сравнивать. В-пятых, зачем вам это? (какая цель у вас для корпуса, который имеет значение)   -  person Yakk - Adam Nevraumont    schedule 30.09.2014
comment
мч, да но к сожалению разницы не вижу. ajclinto, я не думал об этом, наверное, это тоже так. Yakk, мой профессор дал нам задание реализовать различные алгоритмы выпуклой оболочки, я реализовал их все, но понятия не имею, как подойти к ошибкам с плавающей запятой.   -  person jsguy    schedule 30.09.2014
comment
Это не общий ответ, но вы можете использовать два целых числа для точного хранения рационального числа и выполнять всю математику в рациональном пространстве. Это не будет работать, скажем, для sqrt или sin, но будет точным для любого многочлена с рациональными коэффициентами, если независимая переменная также рациональна.   -  person triple_r    schedule 30.09.2014


Ответы (2)


Вы правы в том, что все точки должны быть на выпуклой оболочке, если вы действительно используете точки формы (x, x^2). Однако три точки могут быть коллинеарными. Если вы меняете их или делаете что-то еще странное, это вылетает из окна.

Если вам удастся выбрать свои 100 000 баллов, я бы предложил использовать целые числа в [-50000,49999]. Ваша функция ccw будет вычислять left и right как целые числа, меньшие по модулю, чем 2.5e14 ‹ 2^53, поэтому округления не произойдет.

Ваша сортировка на основе координат будет работать правильно независимо от ввода.

Для общих входных данных следующий предикат ccw содержит ошибку:

inline int ccw(point *p1, point *p2, point *p3) {
  double left  = (p1->x - p3->x)*(p2->y - p3->y);
  double right = (p1->y - p3->y)*(p2->x - p3->x);
  double res = left - right;
  return res > 0;
}

Округление может быть как при вычитании, так и при умножении. Если все ваши точки лежат в ограничивающей рамке H*W, разности координат x будут вычисляться с абсолютной ошибкой около H*eps/2, а разности координат y будут вычисляться с абсолютной ошибкой около W*eps/ 2. Таким образом, продукты будут вычисляться с абсолютной ошибкой около H*W*eps/2. Если fabs(left - right) < 3*H*W*eps/2, нужно точнее оценить left и right. eps здесь 2-52.

Я бы, вероятно, рекомендовал просто использовать MPFR, если сравнение double ничего вам не говорит. Однако вы можете обойтись без него. Трюк с суммированием Кэхана поможет вам получить младшие биты из разностей, а прием 227+1 поможет вам точно вычислить произведения.

person tmyklebu    schedule 30.09.2014

Очень часто в математике с плавающей запятой вам нужно ввести понятие «допуск», иногда обозначаемое как эпсилон. В вашем случае вы можете сделать свою функцию ccw() трехзначной: true/false/indeterminate. Затем, когда вы пытаетесь выяснить, может ли новая точка быть частью выпуклой оболочки, вы спрашиваете: «Это ccw = истина или неопределенность», и в любом случае вы принимаете точку. Неопределенный результат будет иметь место, когда наклон слишком близок к прямой линии, чтобы ее можно было определить.

person John Zwinck    schedule 30.09.2014