Приближенное сравнение чисел в Haskell

[astudent]

«Первое, что приходит в голову» — сравнить подкоренные выражения, а не сами корни.

[Maverik]

Я просто не хотел утяжелять статью ненужными, на мой взгляд, деталями. В оригинале инвариант выглядел так:

prop_mul_distance s p1 p2 =
   distance (s *. p1) (s *. p2) == (abs s) * (distance p1 p2)

Здесь масштабируется не длина вектора, а положение двух точек на плоскости. Функция distance вычисляет расстояние между этими точками, а оператор (*.) умножения на целочисленный скаляр представляет масштабирование координат.

Вот при проверке этого инварианта QuickCheck и обнаружил нарушение. И как теперь, не вскрывая «черного ящика», обеспечить необходимую проверку?

[Komzpa]

Потому в сравнениях для векторов зачастую используют именно квадрат расстояния.

[nerudo]

[Недоверчиво] Точно Haskell виноват?

[Maverik]

Нет, конечно же. Проблема неточности в вычислении глубже, это можно проверить с помощью аналогичного примера на C. Просто в Haskell есть возможность красиво обойти эту проблему. Ну и то, что QuickCheck позволяет найти совершенно неожиданные ситуации, когда что-то работает не так, как ожидалось. Например, (abs s) в инварианте тоже ведь вначале не было.

[eyeless_watcher]

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

Так почему бы не сделать его таковым? Сравнивать, например, вот так:

abs(a-b) < min(abs(a), abs(b)) * 1e-10

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

[Maverik]

Если я правильно понял, это оценка по относительной погрешности? В принципе, неплохая идея. Как-нибудь попробую…

[ZumZoom]

А вы можете заодно, если не сложно, описать пример задачи, в которой необходимо обеспечить совпадение первых десяти значащих разрядов, а не фиксированную точность в шесть (к примеру) знаков после запятой?

[eyeless_watcher]

Она описана в топике последним пунктом. Вычислительные погрешности при использовании IEEE 754 именно такие, пропорциональные масштабу аргументов.

[ZumZoom]

Понял что вы имели ввиду. Кажется, тогда стоит уточнить, в связи с какими соображениями машинный эпсилон изменился на 1e-10. Эта оценка ведь не учитывает вычисления, которые производились с a и b. В общем случае она же не выполняется. Погрешность вполне может перевалить за 1e-10.

[eyeless_watcher]

Изменился от балды. Действительно, не учитывает. Верно, может перевалить. Слово «например» там не просто так :)

[Bodigrim]

Если окажется, что для целей автора достаточно оценивать относительную погрешность, то я бы предложил сделать newtype Circa = Circa Double и переопределить для типа Circa инстанс == необходимым образом. ИМХО это более haskell-way.

[tenzink]

Вопрос не новый. Рекомендую посмотреть в сторону сравнения использующего Units in the last place: здесь. Например, такое сравнение используется в google test framework для сравнения float/double: ASSERT_FLOAT_EQ(expected, actual); означает, что разница <= 4ULP

[cypok]

Нехорошо так делать. Погрешность, о которой вы говорите, на самом деле относительная, а не абсолютная. Будет ли ваш код работать, если x и y будут порядка 1e20?

Лучше будет так:

circaEq t x y = abs (x - y) < t * (x + y)

Но и это не будет работать, если x и y очень близки к нулю, функция усложнится еще больше.

Рекомендую всем, кому это не очевидно, ознакомиться как минимум с The Floating-point Guide, а еще лучше с Comparing Floating Point Numbers, 2012 Edition.

[PsyHaSTe]

Подобно примеру выше, я лично в таком случае пользуюсь просто сравнением относительной разницы, на псевдокоде

circaEq x y precision = 
{
   let mn, mx = min(|x|, |y|), max(|x|, |y|)
   let result = x*y > 0 && max/min - 1 < precision
   result
}

Есть есть фатальный недостаток тут — не стесняйтесь показать.