Пытаясь сделать вариант hypot(x, y), вычисляющий правильно округлённый результат для всех пар чисел x и y, я обнаружил, что используемая для проверки hypot_v3 может быть вовсе не так хороша, как я полагал.

В первой заметке я написал, что она обеспечивает погрешность <= 0.5 ULP, что неточно, так как происходит два округления: при вычислении корня в формате double и затем при конвертировании double во float. Технически погрешность <= 0.5 * (1 + DBL_EPSILON) ULP, эта маленькая добавка (DBL_EPSILON == 2.22045e-16) теоретически может сказаться в редких случаях.

Приведу некоторые очевидные свойства hypot(x, y):

1. Коммутативность: hypot(x, y) == hypot(y, x).
2. Чётность: hypot(-x, y) == hypot(x, -y) == hypot(x, y).
3. Округление: если |x| + |y| == max(|x|, |y|), то hypot(x, y) == max(|x|, |y|).

Данные свойства позволяют оценить количество "разных" пар x, y, дающих нетривиальный результат.

Всего имеется неотрицательных чисел float: 255·2²³ = 2139095040 (включая +бесконечность). Из них субнормальных 2²³-1 = 8388607 (с нулём), нормальных — 2130706433. Для субнормальных имеем 8388607·4194303 = 35184359505921 комбинаций. Для нормальных — 2130706433·2²⁴ = 35747322059030528 комбинаций. Всего 35782506418536449 ≈ 3.58e+16 комбинаций. Если считать, что DBL_EPSILON является хорошей оценкой вероятности ошибки, то умножив её на число комбинаций получим округлённо 8 — матожидание числа неверно округлённых пар. Считанные единицы, на которые можно надеяться натолкнуться только при переборе хотя бы одной восьмой всех возможных сочетаний значений x и y…

Увы, на практике всё несколько хуже. Для нормальных чисел. И много-много хуже для субнормальных! Действительно, оценка погрешности предполагает нормальную форму числа, при которой ULP имеет фиксированный (в относительной погрешности) вес (2^-23). В субнормальных числах последняя цифра имеет больший вес вплоть до 100% (когда это единственная единица — в наименьшем представимом положительном числе). Исходя из этого наблюдения, можно попробовать добавить в hypot_v3 коррекцию порядка:

inline float hypot_v3a(float x, float y)
{
  const double x_ = x, y_ = y,
    s = x_ * x_ + y_ * y_;

  if (s >= 0x1p-251)
    return (float)std::sqrt(s);

  return 0x1p-126f * (float)std::sqrt(0x1p+252 * s);
}

Эта небольшая поправка вроде бы вообще не должна влиять на цифры множителя, однако на практике резко уменьшает число дефектных результатов в субнормальной области (на момент написания я не наткнулся ни на одну подозрительную пару из субнормальных чисел). Среди нормальных подозрение вызывают, например, пары (0.0003162f, 1.661635309e-7f), (0.01f, 0.0001590774482f) и (1e+15f, 4.605317338e+15f).

Итак, я заменил проверочную hypot_v3 на hypot_v3a. Мы знаем, что результат sqrt(x*x + y*y) после коррекции порядка отличается от правильно округлённой гипотенузы не более, чем на вес последней цифры. Здесь надо сделать замечание по поводу того, что именно я под этим подразумеваю. Представление чисел в формате IEEE-754 обладает удобным свойством: их побитовым представлением можно оперировать как целым числом. В случае float 22 младших бита есть двоичные цифры множителя после запятой. Перед запятой подразумевается 1, если число нормальное, и 0, если субнормальное. За 22 битами цифр следует 8 бит порядка. Значение порядка 0 соответствует нулю и субнормальным числам. Значения 1—254 соответствуют нормальным числам (множители от 2^-126 до 2¹²⁷). Порядок 255 соответствует бесконечности (нулевые младшие 22 бита) и нечислам (ненулевые).

В моём случае получается, что модуль разности результата hypot_3a и результата вычисления корня одинарной точности, интерпретируемых как целые числа в соответствии с их двоичным представлением, не превосходит единицы. А значит, получив неточный результат z, я могу рассмотреть два соседних числа, добавив или вычтя 1 из двоичного представления z как из целого числа (то же, что std::nextafter(z, INFINITY), std::nextafter(z, 0.f), но стандартная реализация nextafter тоже не блещет скоростью — она предполагает обработку различных случаев, которых у меня не может быть). Итак, если удастся оценить погрешность результата, то помимо вычисленного значения должно быть достаточно проверить следующее (добавить единицу как к целому) и предыдущее (вычесть единицу как из целого) числа — какое-то из них может оказаться правильным ответом.

В начале, впрочем, я приведу улучшенную версию hypot_v2b, которая избегает работы с субнормальными числами.

inline float hypot_v2d(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_ <= 0x1p-126f)
  {
    invmul = 0x1p-126f;  // 0x1p-23 <= max_ <= 1
    mul = 0x1p+126f;     // 0x1p-23 <= min_ *
  }
  else if (max_ <= 0x1p-63f)
  {
    invmul = 0x1p-87f;   // 0x1p-39 <= max_ <= 0x1p+24
    mul = 0x1p+87f;      // 0x1p-63 <= min_ *
  }
  else if (max_ <= 1.f)
  {
    invmul = 0x1p-24f;   // 0x1p-39 <= max_ <= 0x1p+24
    mul = 0x1p+24f;      // 0x1p-63 <= min_ *
  }
  else if (max_ <= 0x1p+63f)
  {
    invmul = 0x1p+39f;   // 0x1p-39 <= max_ <= 0x1p+24
    mul = 0x1p-39f;      // 0x1p-63 <= min_ *
  }
  else
  {
    invmul = 0x1p+102f;  // 0x1p-38 <= max_ <  0x1p+26
    mul = 0x1p-102f;     // 0x1p-63 <= min_ *
  }

  max_ *= mul;
  min_ *= mul;

  if (max_ + min_ == max_)
    return invmul * max_;
  // * here.
  return invmul * std::sqrt(max_* max_ + min_* min_);
}

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

Заметим, что дискриминант обращается в нуль при , меньше нуля при и больше нуля при . Итак, является правильно округлённым , если , либо и двоичное представление оканчивается на нуль (округление «серединок» к чётным).

Аналогично, выбрав , получим, что больше нуля при .

Теперь воспользуемся double для того, чтобы написать «reference»-версию с оценкой погрешности на основе hypot_v2d.

namespace float_bits
{
  using ifloat = std::int32_t;
  using ufloat = std::uint32_t;
  static_assert(sizeof(ufloat) == sizeof(float));

  inline ufloat to_bits(const float x)
  {
    return reinterpret_cast<const ufloat&>(x);
  }

  inline float to_float(const ufloat bits)
  {
    return reinterpret_cast<const float&>(bits);
  }

  inline bool of_same_sign(float a, float b)
  {
    return ((to_bits(a) ^ to_bits(b)) >> 31) == 0;
  }

  inline float bits_prev(float x)
  {
    return to_float(to_bits(x) - 1);
  }

  inline float bits_next(float x)
  {
    return to_float(to_bits(x) + 1);
  }

  inline ufloat ulp_distance(float a, float b)
  {
    const auto x = to_bits(a), y = to_bits(b);
    return x < y? y - x: x - y;
  }
}

inline float hypot_v2e(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_ <= 0x1p-126f)
  {
    invmul = 0x1p-126f;  // 0x1p-23 <= max_ <= 1
    mul = 0x1p+126f;     // 0x1p-23 <= min_
  }
  else if (max_ <= 0x1p-63f)
  {
    invmul = 0x1p-87f;   // 0x1p-39 <= max_ <= 0x1p+24
    mul = 0x1p+87f;      // 0x1p-63 <= min_
  }
  else if (max_ <= 1.f)
  {
    invmul = 0x1p-24f;   // 0x1p-39 <= max_ <= 0x1p+24
    mul = 0x1p+24f;      // 0x1p-63 <= min_
  }
  else if (max_ <= 0x1p+63f)
  {
    invmul = 0x1p+39f;   // 0x1p-39 <= max_ <= 0x1p+24
    mul = 0x1p-39f;      // 0x1p-63 <= min_
  }
  else
  {
    invmul = 0x1p+102f;  // 0x1p-38 <= max_ <  0x1p+26
    mul = 0x1p-102f;     // 0x1p-63 <= min_
  }

  max_ *= mul;
  min_ *= mul;

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

  const double a = max_, b = min_, s = a*a + b*b;
  const float
    s_ = (float)s,
    s0 = std::sqrt(s_),
    s1 = float_bits::bits_prev(s0),
    s2 = float_bits::bits_next(s0);

  const double d0 = s0, d1 = s1, d2 = s2,
    q0 = s - d0 * d0,
    q1 = s - d1 * d1,
    q2 = s - d2 * d2,
    aq1 = std::fabs(q1),
    aq2 = std::fabs(q2);

  float res = s0;
  if (aq1 < aq2)
  {
    const auto q1_0 = d1 - d0,
      d = 2. * (q1 + q0) + q1_0 * q1_0;
    if (d < 0. || 
       (d == 0. && (float_bits::to_bits(s1) & 1) == 0))
      res = s1;
  }
  else if (aq2 < aq1)
  {
    const auto q2_0 = d2 - d0,
      d = 2. * (q2 + q0) + q2_0 * q2_0;
    if (d > 0. ||
       (d == 0. && (float_bits::to_bits(s2) & 1) == 0))
      res = s2;
  }

  return invmul * res;
}

Результат hypot_v2e совпадает с результатом hypot_v3a почти везде: проверялись все x из набора { 3.16227766e-4f, 0.f, 1e-40f, 0x1p-126f — 0x1p-127f, 0x1p-126f, 1e-30f, 1e-20f, 1e-15f, 1e-6f, 1e-2f, 1.f, 1e+2f, 1e+6f, 1e+15f, 1e+20f, 1e+30f } для всех неотрицательных представимых чисел y — было найдено всего четыре пары (x, y), для которых результаты не совпали — три из них были приведены выше как «подозрительные».

Я написал аналог hypot_v2e, который оперирует только значениями float. Точное возведение в квадрат выполняется с помощью разрезания множителя числа на две 12-битные половинки (hi и lo). Затем вычисляются слагаемые hi*hi, lo*lo и 2*hi*lo. Результат — суммарно 50 отличающихся от hypot_v3a результатов (для того же перебора x-y), т.е. в среднем 1 ошибка на несколько сот миллионов пар (x, y). Тем не менее, не скажу, что такой результат мне особенно понравился. Хотелось бы почти полного совпадения с hypot_v2e или большей простоты кода…