Category Archives: Programming

Один алгоритм раскраски графа II

Недавно я переписал алгоритм субоптимальной правильной вершинной раскраски графа на Python. Решил выложить его сюда.

Вершина графа отождествляется с её индексом. Граф задаётся индексируемой последовательностью (tuple или list) наборов соседей (массив списков соседей). Раскраска — массив (list в терминах Python) цветов.

Итак, функция раскраски у меня имеет следующий вид:

## Построить правильную вершинную раскраску графа.
## \param graph -- граф в виде массива множеств номеров смежных вершин
## \param forbidden -- массив множеств запрещённых цветов для каждой вершины
## \param color -- массив цветов вершин
## \return массив (новых) цветов вершин
def coloring(graph, forbidden=None, color=None):
    if forbidden is None: forbidden = []
    if color is None: color = []
    # Предусловия.
    assert len(forbidden) == 0 or len(forbidden) == len(graph)
    assert len(color) == 0 or len(color) == len(graph)

    if len(forbidden) == 0:
        forbidden = [frozenset()]*len(graph)
    #...

Передав раскраску color, можно задать желаемые цвета, а алгоритм разрешит конфликты, если они есть. Если color уже содержит правильную раскраску, то она и будет возвращена функцией coloring. Массив forbidden позволяет запретить раскрашивать какие-то вершины в какие-то цвета: forbidden[v] есть множество запрещённых цветов (числовых меток) для вершины с индексом v.

Всё, что описано ниже, определено внутри функции coloring.

Если пользователь не передал некую начальную раскраску color, то построим её в два цвета поиском в глубину.

    # Раскраска в два цвета поиском в глубину.
    def coloring2():
        color2 = [0]*len(graph)
        def dfsColoring2(u, c=1):
            if color2[u] == 0:
                for d in range(c, len(graph)+1):
                    if d not in forbidden[u]:
                        color2[u] = d
                        break
                c = 3 - c
                for v in graph[u]:
                    dfsColoring2(v, c)
    
        for start in range(0, len(graph)):
            dfsColoring2(start)
        return color2

    # Построить начальную раскраску, если требуется.
    if len(color) == 0:
        color = coloring2()

Теперь собственно алгоритм. Он заключается в повторении пары действий: выявлении конфликтов цветов (listConflicts) и разрешении конфликтов цветов (resolveConflicts), пока все конфликты не будут разрешены.

    # Перечислить текущие конфликты цветов.
    def listConflicts():
        conflicts = []
        for u in range(0, len(graph)):
            nc = [color[v] for v in graph[u]]
            conflict_rank = nc.count(color[u])
            if conflict_rank > 0:
                nc = frozenset(nc) | forbidden[u]
                new_color = max(nc) + 1
                for c in range(1, new_color-1):
                    if c not in nc:
                        new_color = c
                        break
                conflict = (conflict_rank, new_color, u)
                conflicts.append(conflict)
        conflicts.sort(
            key=lambda c: (c[0], -c[1], c[2]), reverse=True)
        return conflicts
    
    # Разрешить конфликты цветов (потенциально не все).
    def resolveConflicts(conflicts):
        closed = set()
        for _, new_color, u in conflicts:
            if u in closed:
                continue
            color[u] = new_color
            closed |= set(graph[u])
    
    # Собственно алгоритм: начиная с некоторой раскраски,
    # повторять обнаружение и разрешение конфликтов, 
    # пока список конфликтов не опустеет.
    while True:
        conflicts = listConflicts()
        if len(conflicts) == 0:
            return color
        resolveConflicts(conflicts)

Цвет 0 не используется.

Пример:

g = ((1,2), (0,2), (0,1,3,6), (2,4,5), (3,5,6), (3,4), (2,4))
c = coloring(g)
print(c)
> [1, 2, 3, 2, 1, 3, 2]

Раскрашенный граф:

Реклама

Конкатенация и транспонирование 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;
}

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

Пытаясь сделать вариант 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;

Борьба с переполнением и исчезновением при вычислении произведений

Пусть дан массив чисел с плавающей запятой, требуется найти произведение.

double product(const double x[], size_t n)
{
  double p = 1.;
  for (size_t i = 0; i < n; ++i)
    p *= x[i];
  return p;
}

Какие недостатки у приведённого выше очевидного решения?

Из трёх компиляторов (g++ 7.1, clang 4.0, icc 17.0) векторизацию данного цикла (параметры -O3 -mavx2) выполнил только icc, что лишь подтверждает известный тезис о том, что icc порой действует чересчур агрессивно: операция умножения в плавающей точке не ассоциативна.

Для любителей обобщённого кода можно записать то же самое, но для «любых чисел»:

#include <iterator>
#include <functional>
#include <numeric>

template <class It>
auto product(It from, It to)
{
  return std::accumulate(from, to,
    typename std::iterator_traits<It>::value_type(1),
    std::multiplies<>{});
}

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

Одно из возможных решений — складывать логарифмы, затем сумму возвести в степень:

\prod\limits_i{x_i} = \exp\sum\limits_i{\log{x_i}}.

Данный способ, вероятно, лучше подходит для работы с числами с фиксированной запятой (на практике, естественно, лучше использовать логарифм и возведение в степень по основанию 2).

В случае плавающей запятой на основе средств стандартной библиотеки C доступен способ получше: комбинация frexp и ldexp.

Лучше он как с точки зрения скорости (по сути это просто битовые операции над представлением числа), так и с точки зрения точности (не привносит дополнительной погрешности по сравнению с обычным произведением).

#include <cmath>
#include <iterator>

template <class It, class NT>
NT product(It from, It to, NT first)
{
  using std::frexp;
  using std::ldexp;
  int e = 0;
  NT p = frexp(first, &e);
  for (int e2 = 0; from != to; ++from)
  {
    p = frexp(p * *from, &e2);
    e += e2;
  }
  return ldexp(p, e);
}

template <class It> inline
auto product(It from, It to)
{
  return product(from, to,
    typename std::iterator_traits<It>::value_type(1));
}

Или то же самое в виде функтора:

template <class NT, class ET = long>
class Frexp_ldexp_product
{
  NT p;
  ET e;

public:
  using Number = NT;
  using Exponent = ET;

  explicit Frexp_ldexp_product(Number first = Number(1))
  {
    using std::frexp;
    int first_e = 0;
    p = frexp(first, &first_e);
    e = first_e;
  }

  Number operator()() const
  {
    using std::ldexp;
    return ldexp(p, e);
  }

  void operator()(Number next)
  {
    using std::frexp;
    int next_e = 0;
    p = frexp(p * next, &next_e);
    e += next_e;
  }

  Number significand() const
  {
    return p;
  }

  Exponent exponent() const
  {
    return e;
  }
};

// Пример использования: при обработке диапазона.
template <class It, class NT>
NT product(It from, It to, NT first)
{
  return std::for_each(from, to,
    Frexp_ldexp_product<NT>{first})();
}

Функтор удобнее, если значения в произведении генерируются на ходу, а не лежат в памяти. Кроме того, есть возможность получить отдельно мантиссу ∈ (0.5, 1] и порядок, даже если результирующее произведение не помещается в диапазон используемого типа чисел с плавающей точкой.

Напоследок хочется отметить, что если проблема переполнения/исчезновения снята, то рост накопленной погрешности остановить не так просто. Хотя результат произведения зависит от порядка, оценка относительной погрешности от порядка выполнения операций не зависит:

(x(1+n\varepsilon))(y(1+m\varepsilon)) \rightarrow xy(1+n\varepsilon)(1+m\varepsilon)(1+\varepsilon) \rightarrow xy(1+(n+m+1)\varepsilon).

В каком порядке не перемножай, всё равно получится относительная погрешность ne, где n — число сомножителей, e — половина «эпсилона» (ULP).