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

    Наверное, все знают, что при вычислениях с ограниченной точностью два математически эквивалентных выражения могут оказаться не равны друг другу. Например, следующее очевидное математическое равенство при вычислении в 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 экспортирует не только оператор приближенного сравнения, но и комплект функций, для которых он может стать синонимом. Если для приложения неприемлемо использование пико-точности, оно может импортировать необходимую функцию и определить оператор приближенного сравнения соответствующим образом. Точно так же можно поступить, если кому-то больше нравится другое обозначение для оператора приближенного сравнения.
    Share post
    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More
    Ads

    Comments 15

      +2
      «Первое, что приходит в голову» — сравнить подкоренные выражения, а не сами корни.
        0
        Я просто не хотел утяжелять статью ненужными, на мой взгляд, деталями. В оригинале инвариант выглядел так:

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


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

        Вот при проверке этого инварианта QuickCheck и обнаружил нарушение. И как теперь, не вскрывая «черного ящика», обеспечить необходимую проверку?
          +2
          Потому в сравнениях для векторов зачастую используют именно квадрат расстояния.
        +6
        [Недоверчиво] Точно Haskell виноват?
          0
          Нет, конечно же. Проблема неточности в вычислении глубже, это можно проверить с помощью аналогичного примера на C. Просто в Haskell есть возможность красиво обойти эту проблему. Ну и то, что QuickCheck позволяет найти совершенно неожиданные ситуации, когда что-то работает не так, как ожидалось. Например, (abs s) в инварианте тоже ведь вначале не было.
          +6
          Очевидно, что требуемый уровень точности сильно зависит от масштаба чисел, с которыми необходимо работать.

          Так почему бы не сделать его таковым? Сравнивать, например, вот так:
          abs(a-b) < min(abs(a), abs(b)) * 1e-10
          

          И сразу пропадает необходимость в куче функций, отличающихся масштабом.
            0
            Если я правильно понял, это оценка по относительной погрешности? В принципе, неплохая идея. Как-нибудь попробую…
              0
              А вы можете заодно, если не сложно, описать пример задачи, в которой необходимо обеспечить совпадение первых десяти значащих разрядов, а не фиксированную точность в шесть (к примеру) знаков после запятой?
                0
                Она описана в топике последним пунктом. Вычислительные погрешности при использовании IEEE 754 именно такие, пропорциональные масштабу аргументов.
                  0
                  Понял что вы имели ввиду. Кажется, тогда стоит уточнить, в связи с какими соображениями машинный эпсилон изменился на 1e-10. Эта оценка ведь не учитывает вычисления, которые производились с a и b. В общем случае она же не выполняется. Погрешность вполне может перевалить за 1e-10.
                    +2
                    Изменился от балды. Действительно, не учитывает. Верно, может перевалить. Слово «например» там не просто так :)
                0
                Если окажется, что для целей автора достаточно оценивать относительную погрешность, то я бы предложил сделать newtype Circa = Circa Double и переопределить для типа Circa инстанс == необходимым образом. ИМХО это более haskell-way.
                0
                Вопрос не новый. Рекомендую посмотреть в сторону сравнения использующего Units in the last place: здесь. Например, такое сравнение используется в google test framework для сравнения float/double: ASSERT_FLOAT_EQ(expected, actual); означает, что разница <= 4ULP
                  +1
                  Нехорошо так делать. Погрешность, о которой вы говорите, на самом деле относительная, а не абсолютная. Будет ли ваш код работать, если 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.
                    0
                    Подобно примеру выше, я лично в таком случае пользуюсь просто сравнением относительной разницы, на псевдокоде

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


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

                    Only users with full accounts can post comments. Log in, please.