Относительное сравнение чисел с плавающей точкой

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

Допустим, два числа aи b. До сих пор я следил за тем, чтобы как абсолютные ( abs(a-b) < eps), так и относительные ( abs(a-b)/max(abs(a), abs(b)) < eps) ошибки были меньше, чем eps. Таким образом, он улавливает неточности чисел, даже если числа, скажем, около 1e-20.

Однако сегодня я обнаружил проблему, числовое значение aи аналитическое значение bбыли:

In [47]: a                                                                     
Out[47]: 5.9781943146790832e-322

In [48]: b                                                                     
Out[48]: 6.0276008792632078e-322

In [50]: abs(a-b)                                                              
Out[50]: 4.9406564584124654e-324

In [52]: abs(a-b) / max(a, b)                                                  
Out[52]: 0.0081967213114754103

Таким образом, абсолютная ошибка [50] (очевидно) мала, но относительная ошибка [52] велика. Поэтому я подумал, что у меня есть ошибка в моей программе. Отладив, я понял, что эти числа ненормальны . Таким образом, я написал следующую процедуру для правильного сравнения:

real(dp) elemental function rel_error(a, b) result(r)
real(dp), intent(in) :: a, b
real(dp) :: m, d
d = abs(a-b)
m = max(abs(a), abs(b))
if (d < tiny(1._dp)) then
    r = 0
else
    r = d / m
end if
end function

Где tiny(1._dp)возвращает 2.22507385850720138E-308 на моем компьютере. Теперь все работает, и я просто получаю 0 как относительную ошибку, и все в порядке. В частности, приведенная выше относительная ошибка [52] неверна, просто она вызвана недостаточной точностью денормальных чисел. Правильно ли реализована моя rel_errorфункция? Должен ли я просто проверить, что abs(a-b)он меньше, чем крошечный (= денормальный), и вернуть 0? Или я должен проверить какую-то другую комбинацию, как max(abs(a), abs(b))?

Я просто хотел бы знать, что такое «правильный» способ.

Вы можете напрямую проверить денормализованных чисел , используя isnormal()из math.h(C99 или более поздней версии, POSIX.1 или более поздней версии). В Фортране, если модуль ieee_arithmeticдоступен, вы можете использовать ieee_is_normal(). Чтобы быть более точным в отношении нечеткого равенства, вы должны рассмотреть представление ненормальных чисел с плавающей точкой и решить, что вы имеете в виду, чтобы результаты были достаточно хорошими.

Более того, чтобы полагать, что любой результат верен, вы должны быть уверены, что не потеряли слишком много цифр на промежуточном этапе. Вычисления с ненормированными значениями, как правило, ненадежны, и их следует избегать, если ваш алгоритм внутренне масштабируется. Чтобы убедиться, что ваше внутреннее масштабирование прошло успешно, я рекомендую активировать исключения с плавающей запятой, используя feenableexcept()C99 или ieee_arithmeticмодуль в Fortran.

Хотя ваше приложение может перехватить сигнал, возникающий при исключениях с плавающей запятой, все ядра, которые я пробовал, сбрасывают аппаратный флаг, поэтому fetestexcept()не возвращает полезного результата. При запуске с -fp_trapпрограммами PETSc (по умолчанию) распечатывается трассировка стека при возникновении ошибки с плавающей запятой, но она не идентифицирует ошибочную строку. Если вы работаете в отладчике, он сохраняет флаг оборудования и разбивает оскорбительное выражение. Вы можете проверить точную причину, вызвав fetestexceptотладчик, в котором результат является побитовым ИЛИ следующих флагов (значения могут различаться в зависимости от машины, см. fenv.hЭти значения для x86-64 с glibc).

• FE_INVALID = 0x1
• FE_DIVBYZERO = 0x4
• FE_OVERFLOW = 0x8
• FE_UNDERFLOW = 0x10
• FE_INEXACT = 0x20

У Дональда Кнута есть предложение для алгоритма сравнения с плавающей запятой в томе 2 «Получисленные алгоритмы» «Искусства компьютерного программирования». Это было реализовано в C на Th. Belding (см. Пакет fcmp ) и доступен в GSL .

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

Но, приближаясь к нулю, вычислять относительные ошибки бессмысленно.

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

$y$$x$$|y-x|\le absacc+relacc*\max(|x|,|y|)$

Тогда пользователи вашего кода точно знают, насколько они действительно точны.