Author Archive: evetro

Конкатенация и транспонирование tuple

В недавно вышедшей Boost 1.66 beta 1 присутствует новая библиотека метапрограммирования Mp11. Я начал читать отсюда и обратил внимание на два примера.

Первый — конкатенация произвольного набора кортежей типов. Мне показалось очень сложным представленное решение Эрика Ниблера, поскольку «на глаз» было видно, что стандартными средствами это сделать довольно легко, ЕСЛИ ограничиться поддержкой std::tuple. Естественно, я попробовал написать свою версию, которую и привожу ниже.

Вначале конкатенация двух tuple:

template <class A, class B>
struct Tcat;

template <class... T, class... S>
struct Tcat<std::tuple<T...>, std::tuple<S...>>
{
  using type = std::tuple<T..., S...>;
};

template <class A, class B>
using tcat = typename Tcat<A, B>::type;

Теперь не составит труда определить конкатенацию произвольного количества tuple рекурсивно:

template <class... A> struct Tcats;

template <class... T>
struct Tcats<std::tuple<T...>>
{
  using type = std::tuple<T...>;
};

template <class A, class B, class... Other>
struct Tcats<A, B, Other...>
{
  using type = typename Tcats<tcat<A, B>, Other...>::type;
};

template <class... T>
using tcats = typename Tcats<T...>::type;

Второй — transform бинарной операцией двух кортежей (zip, если в терминах функциональных языков программирования). Конкретный случай для двух кортежей не столь интересен. Легко видеть, что при обобщении на произвольное число кортежей мы получаем ни что иное как транспозицию матрицы (F — применяемая операция):


<T11, T12, T13, ..., T1m>,
<T21, T22, T23, ..., T2m>,
...
<Tn1, Tn2, Tn3, ..., Tnm>
--->
F<T11, T21, ..., Tn1>,
F<T12, T22, ..., Tn2>,
F<T13, T23, ..., Tn3>,
...
F<T1m, T2m, ..., Tnm>

Насколько просто сделать такую штуку на шаблонах C++? Результат, который у меня получился, мне не особенно понравился. Думаю, можно это сделать лучше, однако, не факт, что это «лучше» будет означать «короче».

Итак, метафункция, порождающая из столбца исходной матрицы строку новой матрицы с префиксом F (n-арной операцией):

template <class I,
  template <class...> class F, class... T>
struct Tcol
{
  using type = F<std::tuple_element_t<I::value, T>...>;
};

template <class I,
  template <class...> class F, class... T>
using tcol = typename Tcol<I, F, T...>::type;

Здесь I — тип, определяющий зависимое имя value, раскрываемое в целочисленную константу — индекс столбца.

Теперь надо применить tcol к диапазону индексов столбцов. Так как готовых iota, map и bind у меня нет, я написал рекурсивную реализацию аналога «for» «в лоб»:

template <class S, class T,
  template <class...> class F, class... Other>
struct Tfor;

template <class Int, Int s,
  template <class...> class F, class... Other>
struct Tfor<
  std::integral_constant<Int, s>,
  std::integral_constant<Int, s>,
  F, Other...>
{
  using type = std::tuple<>;
};

template <class Int, Int s, Int t,
  template <class...> class F, class... Other>
struct Tfor<
  std::integral_constant<Int, s>,
  std::integral_constant<Int, t>,
  F, Other...>
{
  static_assert(s < t, "");
  using S  = std::integral_constant<Int, s>;
  using S1 = std::integral_constant<Int, s+1>;
  using T  = std::integral_constant<Int, t>;
  
  using type = tcat<std::tuple<F<S, Other...>>,
    typename Tfor<S1, T, F, Other...>::type>;
};

template <class S, class T,
  template <class...> class F, class... Other>
using tfor = typename Tfor<S, T, F, Other...>::type;

Здесь S и T — типы, задающие начало и конец диапазона (::value — соответствующие целочисленные константы), пак Other — исходная матрица типов.

Теперь записать обобщённый zip не составит труда:

template <
  template <class...> class F,
  class... T>
struct Tzip;

template <
  template <class...> class F>
struct Tzip<F> { using type = std::tuple<>; };

template <
  template <class...> class F,
  class T0, class... T>
struct Tzip<F, T0, T...>
{
  template <class I, class... Y>
  using tcolF = tcol<I, F, Y...>;
  using type = tfor<
    std::integral_constant<std::size_t, 
      0>,
    std::integral_constant<std::size_t,
      std::tuple_size<T0>::value>,
    tcolF, T0, T...>;
};

template <
  template <class...> class F,
  class... T>
using tzip = typename Tzip<F, T...>::type;

Пришлось вынуть первый кортеж (строку матрицы) — T0, чтобы узнать количество столбцов. Представляю, какое веселье будет в выводе компилятора, если скормить последовательность кортежей разных размеров :)

Тестировочный код:

template <class... T>
struct Tmax;

template <class T>
struct Tmax<T> { using type = T; };

template <class T0, class... T>
struct Tmax<T0, T...>
{
  using tail_max = typename Tmax<T...>::type;
  using type = std::conditional_t<
    (T0::value < tail_max::value), tail_max, T0>;
};

template <class... T>
using tmax = typename Tmax<T...>::type;

template <class T>
using Sizeof =
  std::integral_constant<std::size_t, sizeof(T)>;

template <class... T>
using tmax_sizeof = tmax<Sizeof<T>...>;

int main()
{
  using namespace std;
  // Concat.
  {
    using t0 = tuple<int, int*>;
    using t1 = tuple<float>;
    using t2 = tuple<char*, char**, void*, void**>;
    using t3 = tuple<>;
    using ex = tuple<
      int, int*, float, char*, char**, void*, void**>;
    static_assert(is_same<tcats<t0, t1, t2, t3>, ex>::value);
  }
  
  // Generic zip.
  {
    using r0 = tuple<char, short, int, long>;
    using r1 = tuple<unsigned, void*, void**, void***>;
    using exr = tuple<
      tuple<char, unsigned>,
      tuple<short, void*>,
      tuple<int, void**>,
      tuple<long, void***>>;
    //print_t(tzip<tuple, r0, r1>{});
    static_assert(is_same_v<exr, tzip<tuple, r0, r1>>);
    
    using t0 = tuple<int[11], int[12], int[13]>;
    using t1 = tuple<int[21], int[22], int[23]>;
    using t2 = tuple<int[31], int[32], int[33]>;
    using t3 = tuple<int[41], int[42], int[43]>;
    using ext = tuple<
      tuple<int[11], int[21], int[31], int[41]>,
      tuple<int[12], int[22], int[32], int[42]>,
      tuple<int[13], int[23], int[33], int[43]>>;
    static_assert(
      is_same_v<ext, tzip<tuple, t0, t1, t2, t3>>);
    
    using exsz = tuple<
      Sizeof<int[41]>, Sizeof<int[42]>, Sizeof<int[43]>>;
    static_assert(
      is_same_v<exsz, tzip<tmax_sizeof, t0, t1, t2, t3>>);
  }
}
Реклама

Вычисление массы молекулы по формуле вещества с помощью Python

Я привык писать на C++, и та задача, решение которой я собираюсь здесь продемонстрировать, сначала была мной решена на C++. Однако Python позволяет решить её с помощью намного более короткой программы (думаю, ряд других популярных языков в этом плане не хуже, но речь сейчас не о них). Это первая запись в моём блоге с использованием Python. Шапочное знакомство с этим языком у меня состоялось много лет назад, но пользоваться им не доводилось. Недавно я решил использовать его на занятиях, в первую очередь, из-за пакета модулей SciPy. В качестве первой задачки для себя я взял задачу, для которой достаточно стандартной функциональности.

Задача. Дана строка, содержащая формулу химического вещества. Требуется найти массу молекулы этого вещества в атомных единицах массы. Формула может содержать названия элементов вида K (одна заглавная латинская буква) или Mg (заглавная + строчная буква), круглые скобки (химики также используют и квадратные скобки, но упростим задачу), натуральные числа (количества атомов и групп атомов). Пример: «Ca (OH)2» (пробелы должны игнорироваться), масса равна сумме масс Ca и удвоенной суммы масс O и H.

Исходные данные: словарь (или файл) с массами элементов.

Первая идея была парсить формулу «руками» подряд, пропуская пробельные символы, выявляя названия элементов и числа, а скобочные формы обрабатывая рекурсивно. Подобная программа на C++ (с комментариями, циклом чтения в main и функцией чтения файла с массами элементов) вылилась примерно в 200 строк кода.

Вторая идея была использовать регулярное выражение для выявления частей формулы. Стандартная библиотека C++ включает поддержку регулярных выражений.

/// Break-up formula into tokens.
vector<string> tokenize_formula(const string &formula)
{
  // Regex object used for tokenizing formulae.
  static const regex tokens("[0-9]+|[A-Z][a-z]?|\\(|\\)");
  vector<string> result;
  sregex_iterator 
    from(begin(formula), end(formula), tokens), to;
  for_each(from, to, 
    [&](auto &match) { result.push_back(match.str()); });
  return result;
}

Дальше всё примерно так же, как при выделении лексем «руками». Такой вариант стал короче почти вдвое, но обзавёлся и важным недостатком — ошибки игнорировались, например, формула «N aHCO3» проходила, как если бы было написано «NHCO3», что неприятно.

Python позволяет записать то же самое ещё короче, но в нём есть возможность ещё сверх того ужать программу, не прибегая при этом к сторонним библиотекам. Причём, если в формуле есть ошибки, то мы это увидим. Впрочем, для начала я приведу вариант на основе предыдущего регэкспа, но уже с отловом ошибок (однако без чтения файла и main). Предполагается, что массы элементов заданы в словаре elements, скажем, так:

elements = {
    'H' :  1.008,
    'C' : 12.01,
    'N' : 14.007,
    'O' : 16,
    'Na': 22.99,
    'Mg': 24.306,
    'P' : 30.974,
    'S' : 32.068,
    'Cl': 35.45,
    'K' : 39.098,
    'Ca': 40.078
    }

Итак, код парсера, используещего методы findall и split:

import re
chemass_tokens = re.compile(r'[0-9]+|[A-Z][a-z]?|\(|\)|\[|\]')

def chemass(formula):
    """
    Формула -- это строка, составленная из знаков элементов, 
    количеств и скобок вроде 'H2O' или 'Ca (OH)2'. 
    Пробелы и всё прочее игнорируется (о выброшенном выводится сообщение).
    @return: масса
    """
    # Получить список токенов:
    tokens = chemass_tokens.findall(formula)
    m, l = do_chemass(tokens)
    # То, что не попало в токены, можно получить,
    # разбив формулу на части, разделённые токенами:
    unparsed = chemass_tokens.split(formula)
    # split возвращает кучу пустых строк, а также пробелы -- выкинем их.
    unparsed = [s for s in unparsed if s != '' and not str.isspace(s)]
    if len(unparsed) > 0:
        print('Chemass ignored the following formula parts:')
        print(unparsed)
    
    return m

def do_chemass(l):
    """
    Отдельная функция нужна для рекурсивной обработки вложенных скобочных форм.
    l -- список строк, где каждая строка -- 
    либо имя элемента, либо число, либо скобка.
    """
    mass = 0
    cm = 0
    while l:
        if str.isdecimal(l[0]):
            cm *= int(l[0])
        elif l[0] == '(' or l[0] == '[':
            mass += cm
            cm, l = do_chemass(l[1:])
            continue # чтобы не удалить голову списка l.
        elif l[0] == ')' or l[0] == ']':
            del l[0]
            break
        elif l[0] not in elements: # сообщить о неизвестном элементе.
            print('Chemass: '+l[0]+' is unknown element, ignored.')
        else:
            mass += cm
            cm = elements[l[0]]
        
        del l[0] # удалить голову списка l.

    mass += cm
    return mass, l

Ну, а теперь, «коронный номер». В Python есть встроенная функция eval, вычисляющая значение выражения (программы) на Python… Преобразуем-ка формулу в выражение, в духе ‘Ca(OH)2’ → ‘Ca+(O+H)*2’ и вычислим его! Для этого можно использовать метод регэкспа sub (подстановка по образцу). Далее привожу полученный код:

import re
chemass_names = re.compile('([A-Z][a-z]?)')
chemass_open = re.compile('(\()')
chemass_num = re.compile('([0-9]+)')

def chemass(formula):
    """
    Формула -- это строка, составленная из знаков элементов, 
    количеств и скобок вроде 'H2O' или 'Ca (OH)2'.
    @return: масса
    """
    return eval(chemass_formula(formula))

def chemass_formula(formula):
    """
    Функция, преобразующая строковую формулу 
    к выражению на Python, вычисляющему массу.
    """
    return chemass_num.sub(r'*\1',
        chemass_names.sub(r'+elements["\1"]',
        chemass_open.sub(r'+\1', formula)))

Это всего 20 строк. Если в формуле есть ошибки, о них сообщит сам парсер Python’а. Да, eval нельзя использовать в некоторых случаях, например, из соображений безопасности, но если это скрипт для личного пользования «по-быстрому», то это может быть очень удобно.

Addendum. Я использовал предварительную поддержку C++17 (g++ 7.2), и меня неприятно удивило отсутствие многих очевидных функций для работы со string_view. Например, его нельзя передать в конструктор fstream или в функцию stoi. Да, всегда можно явно сконвертировать в объект string, но это противоречит цели использования string_view: избежать создания лишних строк.

Проверка возможности записи через итератор вывода

Периодически возникает желание иметь статическую проверку согласованности типа значения итератора вывода и типа записываемого по итератору значения. К сожалению, в том случае, когда в качестве итератора передаётся настоящий итератор вывода, может оказаться, что iterator_traits<OutIt>::value_type есть void. Таковы, в частности, стандартные итераторы вывода, например, back_insert_iterator и ostream_iterator. Причина такого поведения заключается в том, что данные итераторы допускают только присваивание через разыменование и инкремент (который может ничего не делать):

*out = val;

Тип значения, возвращаемого оператором * обычно не является ссылкой на тип записываемого значения (это может быть прокси-объект, в качестве которого может выступать и сам итератор). Подстановка void в качестве value_type приводит к скорым ошибкам компиляции при попытке некорректно использовать этот тип.

Если требуется только проверить доступность такой операции записи для заданного типа значения, то это следует сделать напрямую — проверить возможность указанного выше присваивания:

#include <type_traits>

/// Check if a value type can be written via an iterator.
template <class Val, class It>
using Is_writable_to = typename std::is_assignable<
    decltype(*std::declval<It>()), Val>::type;
     
/// Check if a value type can be written via an iterator.
template <class Val, class It>
constexpr bool is_writable_to 
    = Is_writable_to<Val, It>::value;

Теперь следующий код выводит «110» (double можно привести к int, а void* — нельзя):

#include <iostream>
#include <iterator>
#include <vector>

int main()
{
  using namespace std;
  vector<int> a;
  auto bi = back_inserter(a);
   
  cout << is_writable_to<int, decltype(bi)>;
  cout << is_writable_to<double, decltype(bi)>;
  cout << is_writable_to<void*, decltype(bi)>;
  return 0;
}

Храните деньги в банке

Храните деньги в банке. А банку — под кроватью. В предыдущей записи я неосмотрительно сравнил слово банка (с вариантом баночка) со словом доска, заявив в нём корень банк-. Ну действительно, вроде в русском нет подходящего слова без суффикса (**бана?), да и само слово похоже на недавнее заимствование… Естественно, возникает вопрос, родственны ли слова банк (очевидное позднее заимствование, вероятно, из немецкого Bank) и банка?

Значение «банк» возникло из итальянского banco (видимо, через значения «скамья», «прилавок»), заимствованного из германских языков. Для германского корня Wiktionary.org приводит два слова: *bankiz (откуда англ. bench) и *bankô (откуда англ. bank «берег реки», «отмель»). К модели «изгиб» → «возвышение; вздутие», «возвышение» → «скамья; берег, холм» вполне подходило бы дополнение вида «вздутие» → «медицинская банка» → «банка», в предположении, что медицинские банки пришли в Россию из немецкоговорящих краёв. Но увы в германских языках производные bank- для этой цели, похоже, не используются, а медицинские банки известны весьма продолжительное время.

Ну ладно, посмотрим, что, например, написано в Wiktionary о русском слове: «From Proto-Slavic *banъka, from *banja +‎ *-ъka.» Вот это поворот! Банка — уменьшительное от баня. Но ведь в русском и так есть такое слово: собственно банька. Откуда мог взяться вариант с твёрдым -н-? Сомнительно, чтобы вариант *банъка где-то был зафиксирован непосредственно, но если посмотреть ЭССЯ (т.1, с.152), то, оказывается, в ряде славянских языков есть вариант с мягким -н- и значением, идентичным или близким к «банка». Более того, если посмотреть список слов из предыдущей статьи *ban’a (баня), то там можно найти аналогичное значение («сосуд для воды» или даже «все пузатое») без уменьшительного суффикса. Твёрдость -н- не объясняется, но я могу предложить следующие возможные объяснения: а) в данном значение слово могло быть заимствовано в русский из другого славянского языка, где уже произносилось «твёрдо»; б) когда-то был (но куда-то сплыл) также вариант слова баня с твёрдым исходом; в) повлияло другое слово, например, жбан. Последний вариант можно «вывернуть наизнанку»: чьбанъжбанкабанка под влиянием слова баня, тем более, в ЭССЯ в числе производных праславянского *čьbanъ дано такое слово как слвц. dbanka.

Итак, с этимологической позиции слово банка в значении «сосуд для воды» имеет корень бан-, и я был не прав, поставив его в один ряд со словом доска. (Однако, в русском wiktionary.org указан корень банк-.) Впрочем, следует отметить существование позднего заимствования германского слова «банка» (в значении «скамья» и др.), в котором должен быть корень банк-, но склоняется это слово тоже так, как если бы в нём был суффикс -к-.

Доска, desk, discus

Disclaimer: автор дилетант в лингвистике и просит не судить строго.

Как-то довелось мне встретить рекомендацию писать досочка вместо досточка в качестве уменьшительной формы слова доска. Забавно, мне бы и в голову не пришло писать досточка, но правда ли досочка лучше? Вариант досочка предполагает, что в слове доска корень дос и суффикс к, с чем сложно увязать наличие слов дощатый и дощечка. Ага, вопрос «что правильно досточка или досочка?» некорректен! Правильно — дощечка! Здесь мы имеем чередование в корне: доск/дощ (щ- — результат палатализации ск-). Но точно такое же чередование возможно для корней, оканчивающихся на -ст, что объясняет возникновение варианта досточка то ли из дощатый, то ли из дощечка (депалатализация при замене «смягчающей» -е- на «несмягчающую» -о-). Есть закрепившийся в литературном языке результат аналогичного перехода: пуститьпущатьпускать. Впрочем, есть и примеры закреплённого переосмысления части корня в качестве суффикса (в пользу досочки): банка (корень банк-, upd: см. также) — баночка, а также фляжка (корень фляжк-, ср. англ. flask, нем. Flasche) — фляга (ещё и депалатализация -ж-, почему там вообще -ж-, а не -ш-?). В общем-то, уже форма мн.ч. р.п. досок предполагает переосмысление -к- как суффикса.

Ну ладно, корень доск-. А какова его этимология? Те, кто учил английский, могли заметить подозрительно похожее слово: desk. Этимологи давно заметили это сходство (как и подобные слова в других германских языках), и общим мнением является заимствование из германского *disk(az) славянского *dъska (и в словаре Фасмера, и в ЭССЯ, т.5, с.184), продолжаемого др.-русск. и ст.-слав. дъска. Германское слово заимствовано из латинского discus, а то — в свою очередь, из др.-греческого. Правда, есть технические сложности: германское слово мужского рода и должно было в славянском дать **dьskъ (также другая гласная в корне), чего мы не наблюдаем.

Впрочем, история этого корня в германских языках тоже не так проста: др.-англ. disc [diʃ] не является прообразом desk. В английском есть его прямое продолжение: dish. Это первый «продукт» discus в английском. Третий, поздний, продукт: disc или disk (собственно «диск», с орфографией всё никак не могут определиться). А desk — второй продукт, en.wiktionary.org предоставляет его этимологию: «From Middle English deske, desque, from Medieval Latin desca, modified from Old Italian desco, from Latin discus.» По-моему, это прекрасно. Вариант desca («дэска») вполне годится в качестве прообраза слова дъска. Тем более, что мы встречаем его, например, в ряде славянских языков (по сведениям wiktionary и ЭССЯ). Тут бы хорошо знать момент фиксации слова desca, конечно. Но, думаю, вполне вероятно то, что оно могло появиться много раньше VIII в.

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

Пытаясь сделать вариант 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 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-мерному вектору: найти максимум по абсолютной величине, выбрать порядок, вычислить сумму промасштабированных квадратов с компенсацией и потом уже взять корень, и можно надеяться на погрешность в одну-две последних цифры.

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

Небольшое дополнение к предыдущему посту. Ввиду низкого быстродействия frexp/ldexp может быть разумно заменить их на выбор множителя (степени двойки) банальными if-else.

inline float hypot_v2b(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;
    mul = 0x1p+126f;
  }
  else if (max_ <= 0x1p-63f)
  {
    invmul = 0x1p-63f;
    mul = 0x1p+63f;
  }
  else if (max_ <= 1.f)
  {
    invmul = 1.f;
    mul = 1.f;
  }
  else if (max_ <= 0x1p+63f)
  {
    invmul = 0x1p+63f;
    mul = 0x1p-63f;
  }
  else
  {
    invmul = 0x1p+126f;
    mul = 0x1p-126f;
  }

  max_ *= mul;
  min_ *= mul;

  return invmul * std::sqrt(max_* max_ + min_* min_);
}

Данный код предотвращает некорректное переполнение и исчезновение и обеспечивает погрешность 1 ULP. Конечно, гребёнка if-ов сказывается на скорости работы — у меня получилось, что hypot_v2b примерно втрое медленнее hypot_v1 (но всё же в десять раз быстрее hypot_v2a). Эта оценка довольно условна, так как скорость исполнения ветвлений зависит от множества факторов.

Аналогичный способ можно применить для замены frexp/ldexp при вычислении произведений. Например, если произведение получилось меньше по абсолютной величине наименьшего положительного нормального числа, то пересчитать его, предварительно домножив на 0x1p+63f; если произведение стало бесконечным (переполнение) или приблизилось по величине к наибольшему представимому числу, то пересчитать его, домножив на 0x1p-63f.

Эксперименты с 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.

Вычисление чисел Фибоначчи во время компиляции С++14

Шесть лет назад я привёл пример вычисления целочисленных констант в рамках стандарта C++03. Последний стандарт C++14 ввёл новый вид шаблона: шаблоны переменных. Данная конструкция не только позволяет писать то же самое намного более коротко, но и использовать практически произвольные типы данных.

Пример чисел Фибоначчи, вычисляемых во время компиляции в плавающей точке:

template <unsigned N>
const double fib = fib<N-1> + fib<N-2>;

template <>
const double fib<0> = 0.;

template <>
const double fib<1> = 1.;

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

template <unsigned N>
constexpr unsigned long long fib = fib<N-1> + fib<N-2>;

template <>
constexpr unsigned long long fib<0> = 0;

template <>
constexpr unsigned long long fib<1> = 1;