Tag Archives: sqrt

На пути к вычислению правильно округлённого значения гипотенузы

Пытаясь сделать вариант 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·223 = 2139095040 (включая +бесконечность). Из них субнормальных 223-1 = 8388607 (с нулём), нормальных — 2130706433. Для субнормальных имеем 8388607·4194303 = 35184359505921 комбинаций. Для нормальных — 2130706433·224 = 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 до 2127). Порядок 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_);
}

Пусть есть два соседних числа в плавающей точке a и b = a(1-\varepsilon). Пусть истинное значение корня находится между ними в точке z(t) = (1-t)a + t b. Пограничное значение соответствует t = 1/2. Обозначим разности известных квадратов: \delta_1(t) = z(t)^2 - a^2, \delta_2(t) = z(t)^2 - b^2. Дискриминантом назовём величину

d(t) = 2(\delta_1(t) + \delta_2(t)) + (a - b)^2.

Заметим, что дискриминант обращается в нуль при t = 1/2, меньше нуля при t > 1/2 и больше нуля при t < 1/2. Итак, b является правильно округлённым z(t), если d(t) < 0, либо d(t) = 0 и двоичное представление b оканчивается на нуль (округление «серединок» к чётным).

Аналогично, выбрав b = a(1 + \varepsilon), получим, что d(t) больше нуля при t > 1/2.

Теперь воспользуемся 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 или большей простоты кода…

Реклама

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

В стандартной библиотеке математических функций C++ (унаследованной от C) имеется функция hypot(x, y), возвращающая евклидову длину вектора (x, y).

Зачем особая функция, если можно просто определить «в лоб»?

inline float hypot_v1(float x, float y)
{
  return std::sqrt(x * x + y * y);
}

(Здесь я использую float и только float по причине, о которой будет сказано ниже.)

Причин может быть несколько:

  1. Гарантировать отсутствие переполнения (из-за возведения в квадрат большого числа или сложения больших чисел) в случае, когда конечный результат (длина) представим.
  2. Гарантировать отсутствие влияния исчезновения (из-за возведения в квадрат слишком маленького числа) на результат. В частности, в случае, когда один из аргументов x, y есть нуль, а другой — число, hypot(x, y) должен возвращать модуль другого аргумента, даже если это субнормальное число.
  3. Обеспечивать погрешность лучше 1 ULP.

Легко показать, что hypot_v1 не удовлетворяет всем этим пунктам. Экспериментально установлено, что погрешность hypot_v1 на основной части диапазона, где не возникает ситуация вредоносного переполнения или исчезновения составляет 1 ULP (проверка для y из набора { 0.f, 1e-40f, 1e-30f, 1e-20f, 1e-15f, 1e-6f, 1e-2f, 1.f, 1e+2f, 1e+6f, 1e+15f, 1e+20f, 1e+30f } и для всех возможных значений x).

В активе hypot_v1 сравнительно большая скорость вычисления (у меня получилось, что она примерно в 10 раз быстрее, чем стандартная hypot на g++ -O3).

Можно ли получить какой-то иной компромисс в плане скорости и качества? В Википедии можно наблюдать следующий вариант (на момент написания). Ниже дан исправленный код:

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

  const float
    min_ = x < y? x: y,
    max_ = x < y? y: x;

  if (max_ == 0.f || std::isinf(max_))
    return max_;

  const float q = min_ / max_;
  return max_ * std::sqrt(1.f + q*q);
}

Хорош ли этот код? Он страхует нас от патологического переполнения и исчезновения. Что можно сказать о его скорости? В моём эксперименте получилось, что hypot_v2 примерно вдвое медленнее hypot_v1, что можно считать приемлемым. Что можно сказать о погрешности? Экспериментально получено максимальное значение 2 ULP.

В стремлении убрать операции, вносящие погрешность, я написал иной вариант, использующий вынесение степеней двойки:

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

  if (x == 0.f)
    return y;
  if (y == 0.f)
    return x;

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

  int min_e, max_e;
  min_ = std::frexp(min_, &min_e);
  max_ = std::frexp(max_, &max_e);

  min_ *= min_;
  max_ *= max_;

  return std::ldexp(
    std::sqrt(std::ldexp(min_, 2*(min_e - max_e)) + max_),
    max_e);
}

Этот вариант действительно обеспечивает погрешность <= 1 ULP (экспериментально). К сожалению, этот код в моём случае работал втрое медленнее стандартной hypot. (Что может говорить о потенциально низком быстродействии кода предложенного мной для перемножения массивов чисел.)

Теперь следует прояснить вопрос об экспериментальном вычислении погрешности. Так как используются числа float (IEEE-754 binary32, 23 двоичных разряда после запятой), и доступны числа double (IEEE-754 binary64, 52 двоичных разряда после запятой), то нетрудно организовать вычисление ближайшего к истинной величине значения в формате float, попросту выполнив вычисление с точностью double и округлив затем результат до float:

inline float hypot_v3(float x, float y)
{
  const double x_ = x, y_ = y;
  return (float)std::sqrt(x_ * x_ + y_ * y_);
}

Кстати, этот вариант также не страдает от переполнения и исчезновения и вычисляется лишь вдвое медленнее hypot_v1, обеспечивая погрешность <= 0.5 ULP. Он плох тем, что не масштабируется на double (хотя можно задействовать софтовое вычисление с точностью IEEE-754 binary128, заплатив за это резким падением производительности), зато очень удобен для определения качества других вариантов, включая стандартную hypot.

В случае g++ результат hypot везде совпал с результатом hypot_v3 (т.е. точность <= 0.5 ULP стандартной hypot обеспечена). Интересно, что для её вычисления используется древний (1993г.) код от Sun Microsystems, довольно зубодробительного вида.

В принципе, альтернативой hypot и hypot_v1 может стать их гибрид, использующий hypot_v1 там, где она работает хорошо (с погрешностью 1 ULP), а в остальных случаях (предположительно, редких) вызывающий hypot:

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

  if (1e-19f <= x && x <= 1e+19f &&
      1e-19f <= y && y <= 1e+19f)
    return std::sqrt(x*x + y*y);

  return std::hypot(x, y);
}

Такой код способен работать со скоростью, в среднем близкой к скорости hypot_v1.

О реализации целочисленного аналога функции sqrt

Год назад я написал пост по обращению матриц определённого вида. Он возник не на пустом месте. Исходной посылкой послужило желание сконструировать способ качественного определения быстродействия маленьких быстро выполняющихся функций в условиях, «приближенных к реальности». На тот момент я пытался сравнивать быстродействие различных вариантов реализации функции вычисления целочисленного квадратного корня.

При вычислении sqrt(x) в плавающей точке обычным методом является вычисление 1/sqrt(x) методом Ньютона и затем домножение на x (основная хитрость практической реализации данного метода заключается в способе выбора начального приближения). Эффективно применять этот метод, используя только целочисленные операции, вряд ли возможно, поэтому я решил посмотреть, насколько хорошо (по крайней мере, на конкретном процессоре) в этой задаче работает прямой лобовой метод — двоичный поиск.

Но даже двоичный поиск можно реализовать по-разному. И естественно хочется уметь сравнивать быстродействие разных вариантов в разных условиях.

Неожиданным открытием при первом простом тестировании стала зависимость скорости выполнения функции от её положения в asm-файле (для ассемблерных вариантов). Директива align перед циклами, напротив, не продемонстрировала устойчивого эффекта. Наиболее вероятным «виновником» в данном случае является предсказатель ветвлений и особенности его работы для конкретного CPU. Также может сказываться работа L1I и TLB (здесь, в частности, может удачно или неудачно подействовать директива align).

В случае «простого» тестирования вроде многократного вызова одной функции в цикле есть немалый шанс натолкнуться на проявление таких особенностей во всей их красе. Поэтому я решил поступить несколько хитрее: выполнить псевдослучайное перемешивание разных вызовов. В каждом прогоне — разное количество вызовов разных функций. Выполнив столько же прогонов, сколько у нас есть функций, получим вектор суммарных времён и матрицу количеств вызовов. Решив линейное уравнение, получим оценки времени вызова каждой функции. Можно сделать и больше прогонов и получить решение с помощью МНК (здесь не реализовано). В итоге я всё же не стал использовать аналитически обращаемую матрицу, а задействовал библиотеку uBLAS из состава Boost. См. также о влиянии фрагментации кучи на результаты тестирования.

Использовалась среда Microsoft Visual C++ 2015, компиляция под x86-64 с O2. Причёсывать код однолетней давности я уже поленился — там есть, что выкинуть, но если бы я этим занялся, то… В общем, год уже и так прошёл.

Замечание. В MSVC все исходники .c, .cpp, .asm должны иметь разные имена, так как при компиляции каждого из них порождается одноимённый .obj файл, которые сваливаются в одну общую папку и могут затирать друг друга с печальным результатом во время компоновки.

Итак, исследуемые варианты:

  1. cpp1, asm1 — цикл + ветвление;
  2. cpp2, asm2 — цикл + cmov;
  3. cpp3, asm3 — цикл + ветвление + ilog2 (целочисленный логарифм позволяет выбрать хорошее начальное приближение, для его реализации на x86 используется команда bsr или lzcnt);
  4. cpp4, asm4 — цикл + cmov + ilog2;
  5. asm5 — развёрнутый вариант asm3;
  6. asm6 — развёрнутый вариант asm4 (computed goto);
  7. fpu — использует команды SSE для вычисления isqrt (стандартный метод);
  8. null — ничего не делает, используется для оценки времени вызова.

Итак, результаты тестирования представлены на графике ниже. Числа получены следующим образом: в каждом прогоне использовалось 100’000 псевдослучайных чисел-аргументов, равномерно распределённых от нуля до числа, указанного по оси 0x (1k — 1024, 32k — 32’768, 1M — 220, 1G — 230), выполнялось 1100 повторений, из которых верхние и нижние 5% полученных оценок времён отбрасывались, а от оставшихся 1000 бралось арифметическое среднее. Всего выполнялось по 10 таких прогонов, из результатов которых для каждой функции выбиралось минимальное значение и затем вычиталось соответствующее минимальное значение функции null. В конце все результаты нормировались по функции fpu, так что по оси 0y графиков мы видим, во сколько раз все варианты медленнее вычисления isqrt на fpu.

Процессор, на котором производился тест — AMD K10. Итак, мы видим, что «бодаться» с FPU на «больших» процессорах (пусть даже старых) бессмысленно. Пиррова победа в режиме 1k у asm4/cpp4. Применение ilog2 в комбинации с cmov даёт весьма неплохой эффект, который может оказаться ещё выше на простых процессорах без хорошего предсказателя переходов. А вот разворачивание цикла здесь не помогает (впрочем, опять же на «простом» процессоре ситуация может быть другой).