Эксперименты с hypot III

В посте Эксперименты с hypot я допустил существенную неточность: указал ссылку на реализацию std::hypot (double) без проверки реализации std::hypotf (float), в то время как для сравнения использовалась как раз версия одинарной точности. Итак, hypotf.

Оказывается, это, по сути, та же hypot_v3! Я было думал, что hypotf использует ту же методику, что и hypot. В этой мысли меня поддерживало и её низкое быстродействие (в 5 раз медленнее hypot_v3). Видимо, дополнительные проверки оказались столь дороги. Таким образом, с одной стороны, очевидно, что результаты hypot_v3 и hypotf (для чисел) совпадут. С другой стороны, разработчики руководствовались теми же соображениями, что и я. Функция hypot_v3 даёт корректно округлённый результат в одинарной точности, и её можно использовать для экспериментальной проверки прочих вариантов.

Вариант в двойной точности не мог быть создан с помощью того же приёма, так как четверная точность не поддерживается аппаратно на большинстве процессоров, а программная реализация очень медленная и требует подключения специальной библиотеки (хотя она и есть в комплекте gcc). И вот тут возник вопрос: float vs double даёт возможность проверить качество того или иного алгоритма, а если мы его применяем только в double, то мы, по сути, ограничены теоретическим анализом.

Я, наконец, прочёл комментарий в e_hypot.c, который, вроде бы (с кодом не сравнивал), описывает применённый там метод. А метод основан на соображении следующего вида:

«If (assume round-to-nearest) z=x*x+y*y has error less than sqrt(2)/2 ulp, than sqrt(z) has error less than 1 ulp (exercise).»

Т.е. если вычислить правильно округлённую сумму квадратов, то и корень будет правильно округлён (докажите, ага). Конечно, чисто формально можно разложить в ряд:

\sqrt{z+\varepsilon}-\sqrt{z} = \dfrac{\varepsilon}{2\sqrt{z}} + O(\varepsilon^2).

Если привести абсолютную погрешность к относительной, то получим

\dfrac{\sqrt{z(1+\varepsilon)} - \sqrt{z}}{\sqrt{z}} = \sqrt{1+\varepsilon}-1 = \dfrac{\varepsilon}{2} + O(\varepsilon^2).

В первом приближении извлечение корня уполовинивает относительную погрешность! С такой точки зрения, у нас должен получаться правильно округлённый результат даже в том случае, если последняя цифра z неверна. Впрочем, следует учесть погрешность вычисления самого корня. Если считать, что эпсилон = 0.5 ULP, то у нас появляется дополнительная относительная погрешность \varepsilon:

\left(1 + \dfrac{\varepsilon}{2} + O(\varepsilon^2)\right)(1 + \varepsilon) = 1 + \dfrac{3}{2}\varepsilon + O(\varepsilon^2).

Т.е. всё же ошибка в 0.75 ULP. Я выше написал «правильно округлён», но это, увы, не то же самое, что «ошибка меньше 1 ULP». Ошибку меньше 1 ULP мы получили, но для «правильно округлён» требуется ошибка не более 0.5 ULP, чего мы не достигаем. Эксперимент показывает, что, действительно, вычисление вида

float hypot_v2c(float x, float y)
{
  x = std::fabs(x);
  y = std::fabs(y);

  float
    min_ = x < y? x: y,
    max_ = x < y? y: x,
    mul, invmul;

  if (max_ + min_ == max_)
    return max_;

  // max_ != 0
  // min_ >= 0x1p-24 * max_

  if (max_ <= 0x1p-126f)
  {
    invmul = 0x1p-126f;  // 0x1p-23 <= max_ <= 1
    mul = 0x1p+126f;     // 0x1p-23 <= min_
  }
  else if (max_ <= 0x1p-63f)
  {
    invmul = 0x1p-63f;   // 0x1p-63 <= max_ <= 1
    mul = 0x1p+63f;      // 0x1p-87 <= min_
  }
  else if (max_ <= 1.f)
  {
    invmul = 1.f;        // 0x1p-63 <= max_ <= 1
    mul = 1.f;           // 0x1p-87 <= min_
  }
  else if (max_ <= 0x1p+63f)
  {
    invmul = 0x1p+63f;   // 0x1p-63 <= max_ <= 1
    mul = 0x1p-63f;      // 0x1p-87 <= min_
  }
  else
  {
    invmul = 0x1p+127f;  // 0x1p-64 <= max_ <  2
    mul = 0x1p-127f;     // 0x1p-88 <= min_
  }

  max_ *= mul;
  min_ *= mul;

  const double a = max_, b = min_, s = a*a + b*b;
  return invmul * std::sqrt((float)s);
}

даёт ошибку в последней двоичной цифре, хотя и немного реже, чем вариант hypot_v2b. Следовательно, то же поведение можно ожидать от std::hypot: ошибка в последней цифре. Т.е. она немногим лучше прямого вычисления в духе hypot_v2b. В связи с этим возникает вопрос: можно ли более-менее эффективно вычислить действительно правильно округлённую гипотенузу, не прибегая к арифметике удвоенной разрядности. И ещё: как быть с длиной n-мерного вектора? Метод hypot_v2b можно применить к n-мерному вектору: найти максимум по абсолютной величине, выбрать порядок, вычислить сумму промасштабированных квадратов с компенсацией и потом уже взять корень, и можно надеяться на погрешность в одну-две последних цифры.

Реклама

Добавить комментарий

Заполните поля или щелкните по значку, чтобы оставить свой комментарий:

Логотип WordPress.com

Для комментария используется ваша учётная запись WordPress.com. Выход / Изменить )

Фотография Twitter

Для комментария используется ваша учётная запись Twitter. Выход / Изменить )

Фотография Facebook

Для комментария используется ваша учётная запись Facebook. Выход / Изменить )

Google+ photo

Для комментария используется ваша учётная запись Google+. Выход / Изменить )

Connecting to %s

%d такие блоггеры, как: