Сравнение вещественных чисел

[Blackwanderer]

Широко известно, что сравнивать напрямую числа с плавающей точкой на равенство нельзя:

Код

C++ (Qt)
double a;
double b;
...
if(a == b) // Wrong!!!
...
 

А вот как их сравнивать надо - освещено достаточно плохо. Наиболее распространенный (и дешевый) подход - сравнить разность двух чисел с неким малым значением:

Код

C++ (Qt)
double a;
double b;
...
if(fabs(a - b) < epsilon) // True!!!
...
 

Вот только как выбрать эту самую epsilon? Преимущественное большинство берет значение из головы, в крайнем случае злосчастную DBL_EPSILON. Мне же хотелось бы строго математического обоснования. Кое-что близкое нашел тут и тут, но приведенное обоснование так и не понял (с двоичной арифметикой туговато, надо поднимать). Хотелось бы с вашей помощью разобраться правы или нет авторы материалов по приведенным ссылкам и если правы - то почему. А может вы знаете какой-нибудь еще (более правильный) способ сравнения вещественных чисел?

[alexman]

У qt-ов есть метод, но на чем они основываются хз

Код:

bool qFuzzyCompare(double p1, double p2)

[Blackwanderer]

Код

C++ (Qt)
Q_DECL_CONSTEXPR static inline bool qFuzzyCompare(double p1, double p2)
{
    return (qAbs(p1 - p2) <= 0.000000000001 * qMin(qAbs(p1), qAbs(p2)));
}

Вообще непонятное магическое число.

[Igors]

На эти грабли я наступал много раз Напр нужно вычислить нормаль к треугольнику. На первый взгляд дело нехитрое (треугольник задан 3-мя точками p0, p1, p2)

Код

C++ (Qt)
Point3D nrml = CrossProduct(p1 - p0, p2 - p1);
double len = nrml.length();
if (len < 1.0e-6) ..  
 

Это нормально во многих случаях, напр для рендера: треугольник достаточно мал, значит он НЕ планарен, выкинуть его нафиг - и все дела. Но не всегда. Напр для preview выкидывать его нельзя. Приходится изощряться - попытаться увеличить исходные вектора. Если все равно не помогает - принять нормаль "от фонаря" (напр 0, 0, -1) лучшего нет. Др. словами нет какого-то общего решения, надо смотреть конкретно по задаче.

Конечно есть ф-ции nextafter (чтобы найти следующее значение) но на практике от них мало толку, epsilon почти всегда есть. Существует еще более мерзкая проблема - бой/неустойчивость флота. Ну то если интересно обсудим

[alexman]

Существует еще более мерзкая проблема - бой/неустойчивость флота. Ну то если интересно обсудим

Рассказывай!

[Blackwanderer]

Др. словами нет какого-то общего решения, надо смотреть конкретно по задаче.

Ну вроде как претендует на универсальность следующее:

Код

C++ (Qt)
if (fabs(a - b) <= 16 * DBL_EPSILON * fmax(fabs(a), fabs(b)))
{
    // ...
}
 

[Igors]

Ну вроде как претендует на универсальность следующее:

Ни на что оно не претендует Это довольно прямолинейная попытка вычислить релятивную погрешность. Все зависит от того что Вы собираетесь делать с полученным результатом в задаче. Напр если Вам (хотя бы) потребуется возвести разницу в квадрат - получите проблемы. Более практично установить общий предел точности (напр 1.0e-6) для всех вычислений в задаче - и от него уже плясать.

[Blackwanderer]

Более практично установить общий предел точности (напр 1.0e-6) для всех вычислений в задаче - и от него уже плясать.

Проблема в том, что мне не известен удовлетворительный предел точности. Собственно задача стоит так:

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

Я пробегаюсь по всем треугольникам и нахожу сегменты линии. Неприятности начинаются тогда, когда я их соединяю. А соединяю я их от одного конечного сегмента до другого, присоединяя по одному за раз. Если epsilon слишком большая - то иногда происходит "перескок" через сегмент, т.е. к текущему сегменту подсоединяется не следующий сегмент, а сегмент после следующего. В итоге линия разрывна. Если epsilon слишком маленькая, то иногда два смежных сегмента не считаются смежными. В итоге линия не достраивается.

[Igors]

Ну те "узлы треугольников" называются "вертексы" И нормально надо строить edges (ребра). Т.е. каждый треугольник знает от 0 до 3 смежных треугольников. Тогда ясно какой треугольник брать следующим. Это только на первый взгляд сложновато, но все равно на то и выйдет как только задача хоть немного усложнится. Напр Вы полагаете что у Вас всего одна поверхность, а если их 2 или больше? А если ф-ция пересекает ребро в 2 и более точках? Или просто ф-ция достаточно затратна и лупить по всем треугольникам медленновато. Типичная структура edge

Код

C++ (Qt)
struct CEdge {
 int mVertex[2]  // индексы в контейнере вертексов 
 int mFacet[2];   // индексы треугольников которые шерят это ребро (второй -1 если это край)
};
 
struct CTriangle {
 ...
 int mEdge[3];  // индексы в контейнере CEdge
};
 

[JamS007]

Рассказывай!

+1

Igors, поделитесь, пожалуйста, опытом.

[panAlexey]

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

[kambala]

а че там с qreal? такие же проблемы?

Код

C++ (Qt)
typedef double qreal;

[Igors]

Igors, поделитесь, пожалуйста, опытом.

Ну вот допустим мы удаляем непланарные треугольники (пост #3). И на каждом шаге (или запуске) нам приходят одни и те же исходные треугольники, но повернутые по-разному. Матрица поворота известна, мы можем их развернуть в начальное положение, но при этом данные все же отличаются от исходных (пусть в каком-то дальнем знаке). Этого достаточно чтобы уже другие треугольники удовлетворяли условию планарности. Если есть данные привязанные к ним по индексам - тогда "ой" И можно сколько угодно выбирать epsilon, но всегда найдется такой который "перескочит" через него (и, как правило находится он удивительно скоро )

Поэтому всегда лучше опираться на "порядок построения" чем "склеивать куски" - это ненадежно

[Blackwanderer]

И нормально надо строить edges (ребра). Т.е. каждый треугольник знает от 0 до 3 смежных треугольников. Тогда ясно какой треугольник брать следующим. Это только на первый взгляд сложновато, но все равно на то и выйдет как только задача хоть немного усложнится.

За идею спасибо! В принципе, уже морально подготовился к тому, что придется переделывать все алгоритмы/структуры данных. Вот только процесс это длительный, т.к. озвученная задача - лишь один из модулей программного комплекса. А пока приходится работать с тем что есть. Так что если у кого-то есть мысли/советы/предложения по поводу первоначально поставленного вопроса - милости прошу.

Ну те "узлы треугольников" называются "вертексы"

В терминологии моей предметной области они называются "узлы расчетной области" Я численник.

[Igors]

мысли/советы/предложения по поводу первоначально поставленного вопроса - милости прошу.

epsilon обычно выбирается на основании минимальной длины ребра. Однако читая Ваш пост #7 я не пойму откуда взялась эта проблема? Ну считайте 1 ребро 1 раз - и ничего сравнивать не нужно. Да еще и расчетов в 3 раза меньше  

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