Pull to refresh

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

Reading time 4 min
Views 7.7K
Наверное, все знают, что при вычислениях с ограниченной точностью два математически эквивалентных выражения могут оказаться не равны друг другу. Например, следующее очевидное математическое равенство при вычислении в Haskell неожиданно оказывается ложным:

ghci> 3 * sqrt(24 ^ 2 + 16 ^ 2) == sqrt(72 ^ 2 + 48 ^ 2)
False


Причина такого нарушения в том, что выражения в этом равенстве вычисляются лишь приближенно:

ghci> 3 * sqrt(24 ^ 2 + 16 ^ 2)
86.53323061113574
ghci> sqrt(72 ^ 2 + 48 ^ 2)
86.53323061113575
ghci> sqrt(72 ^ 2 + 48 ^ 2) - 3 * sqrt(24 ^ 2 + 16 ^ 2)
1.4210854715202004e-14


Различие здесь только в последнем (четырнадцатом!) знаке после запятой, но этого уже достаточно, чтобы сравнение оказалось ложным.

Несмотря на то, что эта проблема хорошо известна, программисты уделяют ей мало внимания. Во-первых, считается, что сравнения такого рода возникают только в узкой области численных методов, а во-вторых, что нарушение равенства происходит крайне редко. Как оказалось, и то и другое не совсем верно. Приведенный случай возник, когда мне понадобилось реализовать функцию вычисления длины вектора с целочисленными координатами. При этом для модульного тестирования используются средства пакета QuickCheck, который довольно быстро нашел случай нарушения инварианта масштабирования для длины вектора. Замечу, что это далеко не единственный инвариант, нарушение которого было обнаружено при тестировании.

Возникает вопрос: как проще всего описать проверку приблизительного равенства двух чисел, полученных в результате вычислений с ограниченной точностью? Для решения этой задачи в Haskell достаточно определить еще один оператор сравнения (скажем, ~=), который используется так же, как и обычный оператор равенства. Предлагаю рассмотреть реализацию такого оператора, которую можно оформить в виде достаточно простого модуля Circa.



Первое, что приходит в голову для приближенного сравнения двух чисел — вычислить абсолютное значение разности сравниваемых чисел, и проверить, не превышает ли она некоторого порога:

ghci> abs(sqrt(72 ^ 2 + 48 ^ 2) - 3 * sqrt(24 ^ 2 + 16 ^ 2)) < 1e-12
True


Такое решение, конечно же, вполне работоспособно. Но у него есть два серьезных недостатка. Прежде всего, при чтении этой записи совсем не очевидно, что мы проверяем равенство (пусть даже и приблизительное) двух чисел. Кроме того, нам пришлось воспользоваться «магическим» числом 1e-12, чтобы указать ожидаемую нами неточность сравнения. Конечно же, при небольшом количестве подобных сравнений с этими неудобствами можно было бы смириться. Но когда количество инвариантов измеряется десятками, хотелось бы получить более простой и очевидный способ их записи. Намного лучше выглядит код, в котором двуместный оператор приближенного сравнения используется так же, как и обычный знак равенства, например, вот так:

sqrt(72 ^ 2 + 48 ^ 2) ~= 3 * sqrt(24 ^ 2 + 16 ^ 2)


К сожалению, стандартный синтаксис Haskell не предоставляет такого оператора. Однако сам язык предоставляет замечательную возможность ввести своими силами новый оператор так, как если бы он был частью стандартного синтаксиса.

Для начала определим функцию circaEq, которая сократит запись приближенного сравнения

circaEq :: (Fractional t, Ord t) => t -> t -> t -> Bool
circaEq t x y = abs (x - y) < t


Теперь наш пример становится немного короче, но ненамного понятнее:

ghci> circaEq 1e-12 (sqrt(72 ^ 2 + 48 ^ 2)) (3 * sqrt(24 ^ 2 + 16 ^ 2))
True


Здесь по прежнему много лишней информации: чтобы отделить аргументы, пришлось заключить их в скобки, а самое главное — по прежнему необходимо указывать точность сравнения. Чтобы избавиться от лишнего параметра, воспользуемся трюком, который был применен в модуле Data.Fixed для представления чисел с фиксированной точностью. Определим ряд двуместных функций, каждая из которых будет сравнивать числа с заранее заданной погрешностью. Оказывается, на практике достаточно определить всего семь таких функций для самых типичных значений погрешности:

picoEq :: (Fractional t, Ord t) => t -> t -> Bool
picoEq = circaEq 1e-12
infix 4 `picoEq`

nanoEq :: (Fractional t, Ord t) => t -> t -> Bool
nanoEq = circaEq 1e-9
infix 4 `nanoEq`

microEq :: (Fractional t, Ord t) => t -> t -> Bool
microEq = circaEq 1e-6
infix 4 `microEq`

milliEq :: (Fractional t, Ord t) => t -> t -> Bool
milliEq = circaEq 1e-3
infix 4 `milliEq`

centiEq :: (Fractional t, Ord t) => t -> t -> Bool
centiEq = circaEq 1e-2
infix 4 `centiEq`

deciEq :: (Fractional t, Ord t) => t -> t -> Bool
deciEq = circaEq 1e-1
infix 4 `deciEq`

uniEq :: (Fractional t, Ord t) => t -> t -> Bool
uniEq = circaEq 1
infix 4 `uniEq`


Любая из этих функций может быть использована как двуместный оператор сравнения, например:

ghci> sqrt(72 ^ 2 + 48 ^ 2) `picoEq` 3 * sqrt(24 ^ 2 + 16 ^ 2)
True


Остается всего лишь добавить немного сахара:

(~=) :: (Fractional t, Ord t) => t -> t -> Bool
(~=) = picoEq
infix 4 ~=


и мы получим то, что хотелось:

ghci> sqrt(72 ^ 2 + 48 ^ 2) ~= 3 * sqrt(24 ^ 2 + 16 ^ 2)
True


Ответим еще на один важный вопрос: зачем нам понадобилось несколько разных функций для приближенного сравнения? Неужели сравнения с пико-точностью недостаточно? Оказывается, в самом деле недостаточно. Подходящий контрпример находится с помощью того же пакета QuickCheck:

ghci> sqrt(5588 ^ 2 + 8184 ^ 2) ~= 44 * sqrt(127 ^ 2 + 186 ^ 2)
False
ghci> sqrt(5588 ^ 2 + 8184 ^ 2) `nanoEq` 44 * sqrt(127 ^ 2 + 186 ^ 2)
True


Очевидно, что требуемый уровень точности сильно зависит от масштаба чисел, с которыми необходимо работать. Именно поэтому модуль Circa экспортирует не только оператор приближенного сравнения, но и комплект функций, для которых он может стать синонимом. Если для приложения неприемлемо использование пико-точности, оно может импортировать необходимую функцию и определить оператор приближенного сравнения соответствующим образом. Точно так же можно поступить, если кому-то больше нравится другое обозначение для оператора приближенного сравнения.
Tags:
Hubs:
+10
Comments 15
Comments Comments 15

Articles