Какое сравнение с плавающей точкой является более точным и почему? - PullRequest
Новая
       77

Какое сравнение с плавающей точкой является более точным и почему?

14 голосов
/ 02 марта 2012

Я экспериментирую с различными реализациями метода Ньютона для вычисления квадратных корней. Одним из важных решений является момент прекращения работы алгоритма.

Очевидно, что не стоит использовать абсолютную разницу между y*y и x, где y - текущая оценка квадратного корня из x, потому что для больших значений x это может не можно представить его квадратный корень с достаточной точностью.

Так что я должен использовать относительные критерии. Наивно я бы использовал что-то вроде этого: 

static int sqrt_good_enough(float x, float y) {
  return fabsf(y*y - x) / x < EPS;
}

И это, похоже, работает очень хорошо. Но недавно я начал читать «1013 * элементы стиля программирования» Кернигана и Плаугера , и они дают программу на Фортране для того же алгоритма в главе 1, критерии завершения которого, переведенные в C, будут: 

static int sqrt_good_enough(float x, float y) {
  return fabsf(x/y - y) < EPS * y;
}

Оба математически эквивалентны , но есть ли причина предпочитать одну форму другой?

Ответы [ 2 ]

5 голосов
/ 02 марта 2012

Они все еще не эквивалентны;нижний математически эквивалентен fabsf(y*y - x) / (y*y) < EPS.Проблема, которую я вижу у вас, заключается в том, что если y*y переполнится (вероятно, потому что x - это FLT_MAX и y выбрано неудачно), то прекращение может никогда не произойти.В следующем взаимодействии используются значения типа double. 

>>> import math
>>> x = (2.0 - 2.0 ** -52) * 2.0 ** 1023
>>> y = x / math.sqrt(x)
>>> y * y - x
inf
>>> y == 0.5 * (y + x / y)
True

EDIT: как отмечено в комментарии (теперь удаленном), также полезно разделять операции между итерацией и тестом завершения.

EDIT2: обавероятно, есть проблемы с субнормальным x. профессионалы нормализуют x, чтобы избежать осложнений обеих крайностей.

3 голосов
/ 02 марта 2012

На самом деле они не совсем эквивалентны математически, если только вы не напишите fabsf (y * y - x) / (y * y)

Но я думаю, что ключевым моментом является приведение выражения в соответствие с вашей формулой для вычисления y в итерации Ньютона. Например, если ваша формула y имеет вид y = (y + x / y) / 2, вам следует использовать стиль Кернигана и Плаугера. Если это y = (y * y + x) / (2 * y), вы должны использовать (y * y - x) / (y * y)

Как правило, критерии завершения должны заключаться в том, что abs (y (n + 1) - y (n)) достаточно мала (то есть меньше, чем y (n + 1) * EPS). Вот почему два выражения должны совпадать. Если они не совпадают точно, возможно, что тест завершения решит, что остаток не достаточно мал, в то время как разница в y (n) меньше, чем ошибка с плавающей запятой, из-за различного масштабирования. Результатом будет бесконечный цикл, потому что y (n) перестал изменяться, а критерии завершения никогда не выполняются.

Например, следующий код Matlab - это точно такой же решатель Ньютона, как и в первом примере, но он работает вечно: 

x = 6.800000000000002
yprev = 0
y = 2
while abs(y*y - x) > eps*abs(y*y)
    yprev = y;
    y = 0.5*(y + x/y);
end

В версии на C / C ++ такая же проблема.

...