Category Archives: Lessons

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

Реклама

Пример функциональной декомпозиции

Задача вроде тех, что решают студенты на нашем экзамене после первого семестра.

«Дана прямоугольная матрица. Найти в ней первую строку, не содержащую нулей, или определить, что таких строк нет.»

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

Решение в стиле C подразумевает прототип функции вроде (матрицу можно оформить по-разному, но предпочтём пока простейший вариант):

float* find_row_without_zero
  (float **a, size_t rows, size_t cols);

— вернуть указатель на соответствующую строку, либо nullptr, если такой нет.

Итак, студент видит слово «матрица», из чего делает вывод, что здесь явно двойной цикл, и пишет что-то такое:

float* find_row_without_zero
  (float **a, size_t rows, size_t cols)
{
  for (size_t i = 0; i < rows; ++i)
  {
    for (size_t j = 0; j < cols; ++j)
    {
      if (a[i][j] == 0)
        // ???
      else
        // ???

Почему-то многим очень хочется поставить else-секцию с возвратом результата. Т.е. ошибочно подменить квантор «все» (не нули) на «существует». Под if в этом случае идёт break. Положим, стало ясно, что else здесь никак не нужен. Что же делать? Типичное неказистое решение имеет следующий вид:

// Bad structured programming (using a flag variable).
float* find_row_without_zero_1
  (float **a, size_t rows, size_t cols)
{
  for (size_t i = 0; i < rows; ++i)
  {
    bool found = true;
    for (size_t j = 0; j < cols; ++j)
    {
      if (a[i][j] == 0)
      {
        found = false;
        break;
      }
    }

    if (found)
      return a[i];
  }

  return nullptr;
}

В принципе, кто-то может решиться на использование goto, что даёт избавление от лишней переменной. Но это маловероятно: в примерах goto почти не встречается, а в книгах последних десятилетий его использование традиционно порицается.

// I'm not afraid of GOTO!
float* find_row_without_zero_2
  (float **a, size_t rows, size_t cols)
{
  for (size_t i = 0; i < rows; ++i)
  {
    for (size_t j = 0; j < cols; ++j)
      if (a[i][j] == 0)
        goto next_row;

    return a[i];
  next_row:;
  }

  return nullptr;
}

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

// Good structured programming:
// applying functional decomposition.
float* find_value
  (float *a, size_t size, float val)
{
  for (size_t i = 0; i < size; ++i)
    if (a[i] == val)
      return a + i;
  return nullptr;
}

float* find_row_without_zero_3
  (float **a, size_t rows, size_t cols)
{
  for (size_t i = 0; i < rows; ++i)
    if (!find_value(a[i], cols, 0))
      return a[i];
  return nullptr;
}

Если предполагается использование STL (после первого семестра это ещё не предполагается), то явный цикл можно заменить на вызов std::find, так:

// Ah, std::find does the work in C++.
float* find_row_without_zero_4
  (float **a, size_t rows, size_t cols)
{
  for (size_t i = 0; i < rows; ++i)
    if (std::find(a[i], a[i] + cols, 0.f) == a[i] + cols)
      return a[i];
  return nullptr;
}

или даже так (ещё std::find_if, шаблон по типу данных и значение по умолчанию вместо явного нуля):

// It's not beautiful, I guess...
template <class NT>
NT* find_row_without_zero_5
  (NT **a, size_t rows, size_t cols)
{
  const auto p = std::find_if(a, a + rows,
    [cols](NT *a)
    {
      return std::find(a, a + cols, NT{}) == a + cols;
    });
  return p == a + rows? nullptr: *p;
}

Последний вариант, впрочем, наверно, избыточно «навороченный».

Ну и наконец, можно поразвлекаться в функциональном стиле (требуется C++14), отказавшись от указателей и рассматривая матрицу как диапазон диапазонов. Будем возвращать итератор в стиле STL.

using std::forward;

template <class Item>
auto equals(Item &&item)
{
  return [it = forward<Item>(item)]
  (auto &&other)
  {
    return other == it;
  };
}

template <class Pred>
auto find(Pred &&pred)
{
  return [pr = forward<Pred>(pred)]
  (auto &&rng)
  {
    using std::begin;
    using std::end;

    const auto e = end(rng);
    auto p = begin(rng);
    while (p != e && !pr(*p))
      ++p;

    return p;
  };
}

template <class Pred>
auto none(Pred &&pred)
{
  return [pr = forward<Pred>(pred)]
  (auto &&rng)
  {
    using std::end;
    return find(pr)(rng) == end(rng);
  };
}

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

template <class Rng2>
auto find_row_without_zero_6(Rng2 &&rows)
{
  return find(none(equals(0)))(rows);
}

Забавно, что такая конструкция применима не только, например, к вектору векторов, но и к статическому двумерному массиву:

float m[3][4]
{
  { 1, 1, 1, 0 },
  { 0, 1, 1, 1 },
  { 1, 0, 1, 1 },
};

assert(find_row_without_zero_6(m) == std::end(m));

Обобщённая функция удаления элементов из контейнера

Введение

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

Далее предполагаем следующую «шапку».

#include <algorithm>
#include <iterator>
#include <type_traits>
#include <utility>
using namespace std;

Реализации свободной функции для разных видов контейнеров

Контейнеры-списки реализуют желаемую функцию в виде функции-члена, которой и следует переадресовать вызов.

/// Удаляет элементы, для которых pred возвращает истину,
/// из связных списков (std::list, std::forward_list).
template <class List, class Pred> inline
void list_remove_if(List &l, Pred pred)
{
  l.remove_if(pred);
}

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

/// Удаляет элементы, для которых pred возвращает истину,
/// из линейных контейнеров (std::vector, std::deque).
template <class Vector, class Pred> inline
void vector_remove_if(Vector &v, Pred pred)
{
  v.erase(remove_if(begin(v), end(v), pred), end(v));
}

Ассоциативные контейнеры требуют прямого обхода с удалением конкретных элементов.

/// Удаляет элементы, для которых pred возвращает истину,
/// из ассоциативных контейнеров (разные варианты set и map).
/// В случае map предикат вычисляется для пары ключ-значение.
template <class Assoc, class Pred>
void assoc_remove_if(Assoc &a, Pred pred)
{
  for (auto it = begin(a), it_end = end(a); it != it_end; )
    if (pred(*it))
      it = a.erase(it);
    else
      ++it;
}

 
Автоматический выбор правильного варианта по виду контейнера

Автоматический выбор между представленными тремя вариантами можно организовать через банальную лобовую перегрузку функции remove_if для std::list, std::vector, std::set и т.д. Но мне не нравится этот вариант, так как он не способен автоматически «принять» нестандартные контейнеры, удовлетворяющие интерфейсу того или иного STL-контейнера — для каждого вида контейнера надо будет написать явную перегрузку.

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

Вариант list_remove_if годится тогда, когда контейнер предоставляет метод remove_if (сам контейнер лучше знает как эффективно удалять из него элементы). Вариант vector_remove_if годится для контейнеров с методом удаления диапазона erase(iterator, iterator) и итератором произвольного доступа (формально достаточно итератора однонаправленного последовательного доступа («forward»), но тогда будут подходить и ассоциативные контейнеры). Наконец, вариант assoc_remove_if годится для всех контейнеров с методом iterator erase(iterator).

Далее приводится код метафункций, проверяющих наличие тех или иных операций и таким образом определяющих вид контейнера.

namespace impl
{
  /// Категория итератора для условного типа "не-итератора".
  struct Not_an_iterator_tag {};
  /// Специальный тип, означающий отсутствие типа,
  /// но позволяющий создавать объекты (в отличие от void).
  struct No_type_tag {};
  /// Условный тип "не-итератор", имитирует итератор с
  /// определёнными вложенными типами std::iterator_traits.
  struct False_iterator
    : iterator<Not_an_iterator_tag, No_type_tag> {};
 
  /// Метафункция. Извлекает тип итератора коллекции X 
  /// (использует для этого [std::]begin).
  /// Если begin(X) не определена, возвращает False_iterator.
  template <class X>
  class Collection_iterator
  {
    template <class Y>
    static False_iterator check(unsigned long, Y&&);
 
    template <class Y>
    static auto check(int, Y &&y) -> decltype(begin(y));
 
  public:
    using type = decltype(check(1, declval<X>()));
  };
 
  /// Проверяет, является ли X "коллекцией".
  /// Коллекция -- объект, для которого свободные функции
  /// begin и end возвращают итераторы с категорией,
  /// приводимой к forward_iterator_tag.
  /// Проверяется только [std::]begin.
  /// В качестве коллекции может выступать любой STL контейнер,
  /// объект std::basic_string или статический массив.
  template <class X>
  struct Is_collection
    : is_convertible<typename
        iterator_traits<typename Collection_iterator<X>::type>::iterator_category,
        forward_iterator_tag> {};
 
  ///////////////////////////////////////////////////////////////////////////// 
  /// Проверка, возможен ли вызов assoc_remove_if для объекта Cont.
  template <class Cont>
  class Assoc_remove_if_possible
  {
    template <class C>
    static false_type check(unsigned long, C&&);
 
    template <class C>
    static auto check(int, C &&c) ->
      is_convertible<
        decltype(c.erase(begin(c))),
        decltype(begin(c))>;
 
    static const bool
      erase_defined = decltype(check(1, declval<Cont>()))::value,
      is_collection = Is_collection<Cont>::value;
 
  public:
    using type = integral_constant<bool, erase_defined && is_collection>;
    static const bool value = type::value;
  };
 
  /// Проверка, возможен ли вызов vector_remove_if для объекта Cont.
  template <class Cont>
  class Vector_remove_if_possible
  {
    template <class C>
    static false_type check(unsigned long, C&&);
 
    template <class C>
    static auto check(int, C &&c) -> integral_constant<bool,
         is_convertible<typename
            iterator_traits<typename Collection_iterator<C>::type>::iterator_category,
            random_access_iterator_tag>::value
      && is_convertible<
            decltype(c.erase(begin(c), end(c))),
            decltype(begin(c))>::value>;
 
  public:
    using type = decltype(check(1, declval<Cont>()));
    static const bool value = type::value;
  };
 
  /// Проверка, возможен ли вызов list_remove_if.
  template <class Cont>
  class List_remove_if_possible
  {
    template <class C>
    static false_type check(unsigned long, C&&);
 
    // Универсальный предикат.
    struct Unipred
    {
      template <class T> 
      bool operator()(const T&) { return false; }
    };
 
    template <class C>
    static auto check(int, C &&c) ->
      decltype((void)c.remove_if(Unipred{}), true_type{});
 
  public:
    using type = decltype(check(1, declval<Cont>()));
    static const bool value = type::value;
  };
 
  ///////////////////////////////////////////////////////////////////////////// 
  // Статическая диспетчеризация вызова
  // Теги для статического выбора соответствующего варианта remove_if.
  struct Assoc_remove_if_tag {};
  struct Vector_remove_if_tag {};
  struct List_remove_if_tag {};
 
  /// Выбор подходящего тега remove_if для контейнера.
  template <class Cont>
  class Remove_if_tag
  {
    static const bool
      assoc_possible  = Assoc_remove_if_possible<Cont>::value,
      vector_possible = Vector_remove_if_possible<Cont>::value,
      list_possible   = List_remove_if_possible<Cont>::value;
 
    static_assert(assoc_possible || vector_possible || list_possible,
        "remove_if can not be called for this type");
 
  public:
    using type = conditional_t<list_possible, List_remove_if_tag,
                 conditional_t<vector_possible, Vector_remove_if_tag,
                 conditional_t<assoc_possible, Assoc_remove_if_tag,
                               void>>>;
  };

  // Перенаправление вызова remove_if на варианты для различных типов контейнеров.
  template <class Cont, class Pred> inline
  void remove_if(Cont &c, Pred pred, Assoc_remove_if_tag)
  { assoc_remove_if(c, pred); }
 
  template <class Cont, class Pred> inline
  void remove_if(Cont &c, Pred pred, Vector_remove_if_tag)
  { vector_remove_if(c, pred); }
 
  template <class Cont, class Pred> inline
  void remove_if(Cont &c, Pred pred, List_remove_if_tag)
  { list_remove_if(c, pred); }
}

/// Вызов разрешается статически в подходящий вариант в зависимости от типа контейнера.
template <class Cont, class Pred> inline
void remove_if(Cont &container, Pred pred)
{
  impl::remove_if(container, pred, typename impl::Remove_if_tag<Cont>::type{});
}

Итак, remove_if(контейнер, предикат) выбирает один из трёх способов удаления элементов, удовлетворяющих предикату, в зависимости от операций, предоставляемых контейнером. Приведённый код целиком и небольшой проверочный код: ideone, Bitbucket.

UPD. Увидел, что нечто подобное предлагалось для добавления в стандартную библиотеку. Пока оно представлено в виде расширения в пространстве имён experimental.

О реализации одного алгоритма топологической сортировки

Как-то так получилось, что именно записи на в общем-то побочную для меня тему алгоритмов на графах привлекают основной объём посещений блога «со стороны». Отчасти поэтому, а также потому, что на лекции я предложил было не слишком разумный вариант реализации одного из алгоритмов топологической сортировки вершин ориентированного ациклического графа (DAG — «directed acyclic graph»), я решил написать данный текст и соответствующий код.

На изображении ниже представлен DAG, использованный для проверки реализованной программы.

Пример DAG

Пример DAG

Рассматриваемый алгоритм весьма прост. Его можно описать следующим образом.

Вход: граф G = (V, E).
Выход: последовательность R.
Пусть R = ∅.
Пока |R| < |V|:
    пусть v ∈ V \ R такая что indeg(v) = 0;
    дописать v к R;
    удалить v из G.

Хотя этот алгоритм описан как «деструктивный» по отношению к графу G, на практике совсем не обязательно «разрушать» пусть даже копию некоторого представления G. Его можно переформулировать следующим образом (предполагая представление графа в виде списков соседей).

Создать A = [ indeg(v) | v ∈ V ] — массив степений захода вершин.
Завести очередь Q и поместить туда все v ∈ V: A[v] = 0;
Пока Q ≠ ∅:
    извлечь очередную v из Q;
    поместить v в R;
    для всех w ∈ out(v):
        если A[w] = 1: поместить w в Q;
        если A[w] > 0: A[w] := A[w] – 1.

Как нетрудно видеть, данный алгоритм имеет сложность O(|V| + |E|), при предположении, что затраты на перечисление соседей вершины пропорциональны их количеству. Никакие «сложные» структуры данных не используются.

Ссылка на реализацию на C++.

Метод Эстрина

Недавно мне попался на глаза интересный метод вычисления многочленов, называемый в источнике методом Эстрина (предполагаю, что имеется в виду Джеральд Эстрин, но я глубоко не копал). Презентацию по приведённой ссылке рекомендуется прочитать тем, кого интересует работа с плавающей точкой и вычисление значений функций. Здесь же я рассмотрю лишь небольшой фрагмент.

Итак, пусть дан массив чисел, представляющий набор коэффициентов многочлена по целым степеням x. Я буду обозначать его a, а его размер – n.

Создадим функцию, вычисляющую для заданных набора коэффициентов и значения переменной x значение многочлена P(x). Боюсь, большинство студентов напишет примерно следующее

double poly(double x, double a[], int n)
{
  double s = 0;
  for (int i = 0; i < n; ++i)
    s += a[i] * pow(x, i);
  return s;
}

Выдав такое, например, на собеседовании, кандидат рискует быть отправленным куда подальше.

Запретив вызывать дорогую функцию pow, мы получим, скорее всего, следующий вариант

double poly(double x, double a[], int n)
{
  double s = 0, xi = 1;
  for (int i = 0; i < n; ++i)
  {
    s += a[i] * xi;
    xi *= x;
  }
  return s;
}

Впрочем, может оказаться, что студент учился неплохо и вспомнит тот стандартный метод вычисления многочленов, который мы справедливо ожидали бы увидеть и на собеседовании. Метод Горнера. Данный метод позволяет сэкономить одно умножение на каждой итерации цикла суммирования. Будем считать, что пустой массив нам не передадут: в конце концов, пустой массив коэффициентов означает, что функция, которую мы пытаемся вычислить, не определена. И кстати, заменим-ка int…

double poly(double x, double a[], size_t n)
{
  assert(n != 0);
  double s = a[--n];
  while (n != 0)
    s = s * x + a[--n];
  return s;
}

Движение назад по массиву коэффициентов несколько неприятно, но с этим можно смириться.

Метод Горнера основан на часто встречающейся тернарной операции, называемой умножение-сложение (multiply-add «mad»/«madd», multiply-accumulate «mac»/«macc»). Многие процессоры (CPU и GPU) поддерживают данную операцию на уровне машинных команд (нередко и для целых чисел тоже). Что до плавающей точки, то существует два варианта реализации: с двумя округлениями (последовательное умножение и сложение) и с одним округлением (слитое умножение-сложение, fused multiply-add, «fma»). Второй вариант может быть довольно интересен с точки зрения уменьшения погрешности вычислений и был добавлен в стандарт IEEE-754 2008 года. Процессоры семейства x86 с недавних пор располагают аппаратной реализацией данного действия (наборы команд FMA4, FMA3). Программная реализация возможна, но будет на порядок проигрывать аппаратной с точки зрения быстродействия (которое для машинных команд FMA близко или равно быстродействию обычного умножения). Далее, стандарт языка C принял FMA ещё раньше: ISO C99 math.h включает функции fma, fmaf, fmal, которые могут задействовать соответствующие ресурсы процессора, либо быть основаны на программной реализации. ISO C++11 также предполагает наличие fma в cmath (в MS Visual C++ есть поддержка, начиная с версии 12 (2013)). Аппаратная поддержка fma открывает возможность для эффективных программных реализаций операций деления и взятия квадратного корня с разной степенью допустимой погрешности.

Заметим, что компилятор по умолчанию не имеет права заменять выражения вида (a*b+c) на инструкцию fma или вызов аналогичной функции, так как это влияет на результат (разная погрешность). Поэтому, если мы хотим меньшую погрешность данной операции, то следует явно вызывать функцию fma. Если мы хотели бы меньшую погрешность, но не готовы платить за это резким падением скорости вычислений, то C99 предоставляет макрос FP_FAST_FMA, который должен быть определён, если для вычисления fma используется аппаратная реализация.

Итак, используя fma, получаем вариант

double poly_fma(double x, double a[], size_t n)
{
  assert(n != 0);
  double s = a[--n];
  while (n != 0)
    s = fma(s, x, a[--n]);
  return s;
}

К сожалению, данная реализация далеко не идеальна с точки зрения быстродействия. Современные процессоры обычно включают два и более конвейеров, способных выполнять инструкции, оперирующие значениями в плавающей точке. Однако, чтобы загрузить их работой, выполняемые инструкции должны быть максимально независимы друг от друга. В методе Горнера каждая следующая команда fma использует результат предыдущей (сумму s). Ситуация обостряется при возможности задействовать SIMD-команды, например, Intel SSE или ARM NEON – метод Горнера, увы, сугубо последователен.

Метод Эстрина использует ассоциативность и дистрибутивность операций сложения и умножения, чтобы построить сбалансированное дерево вычисления многочлена. На схеме ниже приведён пример такого дерева для размера восемь (седьмая степень).

Метод Эстрина для многочлена седьмой степени

Это соответствует следующему коду

const auto
  x2 = x * x,
  a10 = fma(a[1], x, a[0]),
  a11 = fma(a[3], x, a[2]),
  a12 = fma(a[5], x, a[4]),
  a13 = fma(a[7], x, a[6]),
  x4 = x2 * x2,
  a20 = fma(a11, x2, a10),
  a21 = fma(a13, x2, a12),
  result = fma(a21, x4, a20);

Параллельное вычисление ветвей дерева позволяет эффективно загрузить конвейеры CPU, либо использовать SIMD-операции. Сложно сказать, насколько эффективна будет, например, реализация данного метода для произвольного массива в рекурсивной форме – накладные расходы могут перевесить. Поэтому при работе с произвольным многочленом представляется разумным явно вычислять лишь нижние уровни дерева, суммируя их результаты тем же методом Горнера, но уже для некоторой степени x (квадрат для нижнего уровня, четвёртая степень для двух нижних уровней и т.д.).

Например, если раскрыть только самый нижний уровень дерева и явно вычислять комбинации вида (a[2*i+1]*x + a[2*i]), то получим примерно следующий код

double poly_estrin2_fma(double x, double a[], size_t n)
{
  assert(n != 0);
  if (n == 1)
    return a[0];

  // один шаг схемы Эстрина + метод Горнера + fma
  const auto x2 = x * x;
  auto s = fma(a[n-1], x, a[n-2]);
  for (--n; n > 2; n -= 2)
  {
    const auto p = fma(a[n-2], x, a[n-3]);
    s = fma(s, x2, p);
  }

  // n==2, если был нечётный размер массива
  return n != 2? s: fma(s, x, a[0]);
}

Но по теоретическим прикидкам получается, что для того чтобы загрузить два конвейера, следует раскрыть хотя бы два нижних уровня дерева. Я выразил это в более обобщённом виде с возможностью выбирать способ вычисления умножения-сложения через параметр шаблона.

struct Multiply_add
{
  template <class Float>
  static Float madd(Float a, Float b, Float c)
  {
    return a * b + c;
  }
};

struct Fused_multiply_add
{
  template <class Float>
  static Float madd(Float a, Float b, Float c)
  {
    using std::fma;
    return fma(a, b, c);
  }
};

template <class Madd = Multiply_add>
class Estrin_4 : Madd
{
  using Madd::madd;
public:

  template <class Float, class Float_it>
  static Float compute(Float x, Float_it a, size_t n)
  {
    assert(n != 0);
    Float s, x2 = x * x;
    // начало считаем по разным формулам в зависимости от остатка
    const size_t remainder = n % 4;
    switch (remainder)
    {
    case 1: s = a[n-1]; break;
    case 2: s = madd(a[n-1], x, a[n-2]); break;
    case 3: s = madd(madd(a[n-1], x, a[n-2]), x, a[n-3]); break;
    default:
      s = madd(madd(a[n-1], x, a[n-2]), x2, madd(a[n-3], x, a[n-4]));
      n -= 4;
    }

    if (n -= remainder)
    {
      const auto x4 = x2 * x2;
      do
        s = madd(s, x4,
              madd(madd(a[n-1], x, a[n-2]), x2,
                madd(a[n-3], x, a[n-4])));
      while (n -= 4);
    }

    return s;
  }
};

// пример использования класса
inline double poly_estrin4_fma(double x, double a[], size_t n)
{
  return Estrin_4<Fused_multiply_add>::compute(x, a, n);
}

Beware! Я провёл несколько тестов данного кода с многочленами разного размера и он дал правильный результат, однако я не даю никаких гарантий относительно его корректности и, тем более, реального быстродействия (это я не проверял).

C++: возможная замена стандартным min & max

Стандартный заголовок algorithm предоставляет определения min, max и minmax (возвращает std::pair). Итак, как выглядит объявление, например, базового варианта std::max?

template <class T>
const T& max(const T& a, const T& b);

Недостатки этого объявления:

  • одинаковый тип a и b — наиболее раздражающий недостаток. Например, что-нибудь безобидное вроде max(vec.size(), 100) уже не компилируется, потому что 100 имеет тип int, а не size_t;
  • только два параметра, поэтому max(a, b, c) придётся оформить как max(a, max(b, c));
  • передача по const-ссылке:
    • для переменных a и b невозможно написать max(a, b) = c;
    • невозможна перемещающая семантика для rvalue: auto v = max(compute_vector1(), compute_vector2()) — всегда копирование;
    • иногда возникает проблема с порождением оптимального кода для значений простых типов.

Что до недостатка 2, то C++2011 предоставляет новые варианты, принимающие initializer_list:

template <class T>
T max(std::initializer_list<T> ilist);

Теоретически такой вариант ликвидирует все недостатки, в том числе, может предоставлять перемещающую семантику, но я не изучал данный вопрос подробно (в msvc2013 выполняется копирование).

Рассмотрим недостаток 1. Простое определение для двух разных типов вызывает затруднение в месте описания возвращаемого типа. C++2014 мог бы позволить написать просто

template <class U, class V> inline
auto max(U u, V v)
{ return u < v? v: u; }

В C++2011 придётся прибегнуть к decltype:

template <class U, class V> inline
auto max(U u, V v)
  -> decltype(u < v? v: u)
{ return u < v? v: u; }

или воспользоваться метафункцией std::common_type из type_traits (также доступна в варианте Boost для ряда старых компиляторов):

template <class U, class V> inline
typename std::common_type<U, V>::type
max(U u, V v)
{ return u < v? v: u; }

Попробуем обобщить данный код, используя универсальные ссылки (std::forward из memory):

template <class U, class V> inline
typename std::common_type<U, V>::type
max(U &&u, V &&v)
{
  return u < v? std::forward<V>(v): std::forward<U>(u);
}

К сожалению, код

int a = 1, b = 2;
max(a, b) = 3;

не компилируется, потому что common_type убирает ссылки (возвращаемый тип данного варианта max будет int, а не int&). Вероятно, простейший вариант исправить этот недостаток — реализовать объединение типов с помощью оператора ?: (сохраняющего ссылочные типы).

template <class U, class V> inline
auto max(U &&u, V &&v)
  -> decltype(u < v? std::forward<V>(v): std::forward<U>(u))
{
  return u < v? std::forward<V>(v): std::forward<U>(u);
}

Обобщим это определение на произвольное число параметров.

template <class U, class V, class ...Tail> inline
auto max(U &&u, V &&v, Tail&&... tail)
  -> decltype(max(std::forward<U>(u),
       max(std::forward<V>(v), std::forward<Tail>(tail)...)))
{
  return max(std::forward<U>(u),
         max(std::forward<V>(v), std::forward<Tail>(tail)...));
}

К сожалению, этот код не компилируется (при попытке вычислить max(1, 2, 3, 4), например) компилятором g++ (хотя компилируется msvc2013). Проблема в рекурсивной зависимости вывода типа — при вычислении max от 4 аргументов вариант для 3-х ещё не сгенерирован. Поэтому желательно отделить вычисление возвращаемого типа в метафункцию.

При использовании имён max/min вместе с объектами стандартных классов возможна «кража функций» (function hijacking) механизмом ADL (а то и срабатывание макросов max и min), поэтому я буду использовать (пусть и непривычные) имена maxi и mini, заведомо не конфликтующие со стандартными определениями. Итак, на примере минимума:

template <class...>
struct commonref_type {};

template <class U>
struct commonref_type<U>
{
  using type = U;
};

template <class U, class V>
struct commonref_type<U, V>
{
private:
  static U make_u();
  static V make_v();
  static bool dummy_condition();

public:
  using type = decltype(dummy_condition()? make_u() : make_v());
};

template <class U, class Head, class ...Tail>
struct commonref_type<U, Head, Tail...>
{
  using type = typename
    commonref_type<U, typename
    commonref_type<Head, Tail...>::type>::type;
};


template <class U, class V> inline
typename commonref_type<U, V>::type
mini(U &&u, V &&v)
{
  return v < u? std::forward<V>(v) : std::forward<U>(u);
}

template <class U, class V, class ...Tail> inline
typename commonref_type<U, V, Tail...>::type
mini(U &&u, V &&v, Tail&&... tail)
{
  return mini(std::forward<U>(u),
         mini(std::forward<V>(v), std::forward<Tail>(tail)...));
}

Формат, позволяющий передать произвольное количество разнотипных параметров, закрывает возможность с помощью механизма перегрузки функций добавить варианты, принимающие произвольный компаратор (такие есть у стандартных min и max). Но ничто не мешает дать другие имена, например, minc и maxc (от comparator). В отличие от стандартных вариантов, здесь будем передавать компаратор первым параметром.

template <class Comp, class U, class V> inline
typename commonref_type<U, V>::type
minc(Comp &&comp, U &&u, V &&v)
{
  return comp(v, u)? std::forward<V>(v) : std::forward<U>(u);
}

template <class Comp, class U, class V, class ...Tail> inline
typename commonref_type<U, V, Tail...>::type
minc(Comp &&comp, U &&u, V &&v, Tail&&... tail)
{
  return minc(comp, std::forward<U>(u),
         minc(comp, std::forward<V>(v), std::forward<Tail>(tail)...));
}

Последнее, с чем осталось разобраться — реализация «одновременного» поиска минимума и максимума. Следуя стандартному определению, будем возвращать пару (результат min, результат max) [std::tuple].

template <class ...Types>
struct minmax_result
{
  using part = typename commonref_type<Types...>::type;
  using type = std::tuple<part, part>;
};

Корректная реализация функций одновременного выбора минимума и максимума требует большого внимания к деталям. Моя версия:

template <class U, class V> inline
typename minmax_result<U, V>::type
minmaxi(U &&u, V &&v)
{
  using RT = typename minmax_result<U, V>::type;
  return v < u?
    RT(std::forward<V>(v), std::forward<U>(u))
  : RT(std::forward<U>(u), std::forward<V>(v));
}

template <class U, class V, class ...Tail> inline
typename minmax_result<U, V, Tail...>::type
minmaxi(U &&u, V &&v, Tail&&... tail)
{
  using RT = typename minmax_result<U, V, Tail...>::type;
  using R = typename commonref_type<V, Tail...>::type;
  auto tail_result =
    minmaxi(std::forward<V>(v), std::forward<Tail>(tail)...);
  auto &p0 = std::get<0>(tail_result),
       &p1 = std::get<1>(tail_result);
  return u < p1?
        (p0 < u?
             RT(std::forward<R>(p0), std::forward<R>(p1))
           : RT(std::forward<U>(u),  std::forward<R>(p1)))
           : RT(std::forward<R>(p0), std::forward<U>(u));
}

Текущий вариант реализации целиком.

Метафункции C++

Общее определение: функции, вычисляемые во время компиляции, в том числе, позволяющие вычислять типы.

Механизм шаблонов C++ является чистым функциональным языком (с громоздким синтаксисом), выражения на котором вычисляются во время компиляции, поэтому основной метод реализации метафункций опирается на механизм шаблонов и (частичной) специализации. В C++2011 появился второй метод: на основе constexpr, здесь он не рассматривается.

Я буду выписывать «определения» в псевдокоде, напоминающем Haskell, а затем приводить аналог на C++.

Связывание имени с константой

let x = y
// x -- новое имя для значения y
// применительно к типам в C++
// C++1998
typedef y x;

// C++2011
using x = y;
Вариант 1

Пусть f принимает тип и возвращает тип.

Определение

f :: type -> type
f x = ...
template <class x> class f { ... };

Применение

f x  // или
f(x) // или
(f x)
f<x> // C++

Функции высших порядков

f :: (type -> type) -> type
f x = ...
template <template <class> x> class f { ... };

Достоинство: краткость синтаксиса применения.

Недостатки:

  • нельзя сопоставить уже существующие типы (например, нельзя выразить определение f int = float);
  • неудобно использовать функции высших порядков.
Вариант 2

Прообраз введён стандартной библиотекой C++, использующей зависимые имена в шаблонах (например, в std::unary_function, классах характеристик std::char_traits, std::iterator_traits, контейнерах).
Популяризован библиотеками Boost.MPL, Boost.TypeTraits (последняя вошла в стандарт C++2011).

Определение

f :: type -> type
f x = expr
template <class x> struct f
{
  typedef expr type;
};

Применение

f x // f(x)
typename f<x>::type // C++

Достоинства:

  • можно привязать любой тип;
  • можно задать мультифункцию (несколько зависимых имён);
  • можно вычислять не только типы, но и значения (традиционно используется вложенное имя value):
factorial :: unsigned -> unsigned
factorial n = n * factorial (n - 1)
factorial 0 = 1
template <unsigned n> struct factorial
{
  enum { value = n * factorial<n - 1>::value };
};

template <> struct factorial<0> { enum { value = 1u }; };

Другой классический пример такого рода: числа Фибоначчи.

  • можно задавать через наследование (используется std::conditional из type_traits):
uint :: unsigned -> type
uint bits = if bits < 9  then uint8_t  else
            if bits < 17 then uint16_t else
            if bits < 33 then uint32_t else
            uint64_t ;
template <unsigned bits> struct uint
  : conditional<(bits < 9),  uint8_t,  typename
    conditional<(bits < 17), uint16_t, typename
    conditional<(bits < 33), uint32_t,
    uint64_t >::type>::type> {};

Недостатки:

  • длинный синтаксис применения (причём typename требуется в контексте определения шаблона при использовании имён, зависящих от параметров шаблона, но запрещено в противном случае);
  • неудобно передавать в функции высшего порядка.

Замечание
С++2011 позволяет сократить синтаксис применения с помощью template-using объявления:

template <unsigned bits> using uint_t = typename uint<bits>::type;
uint_t<19> counter; // OK

Данный подход применяется в type_traits C++2014.

Замыкания и карринг

mul a b = a * b
mul5 = mul 5
thirty = mul5 6

Карринг в C++ напрямую не поддерживается. Его можно имитировать, разделив параметры шаблона.

template <int a>
struct mul
{
  template <int b>
  struct curried { enum { value = a * b }; };
};

using mul5 = mul<5>;
const auto thirty = mul5::curried<6>::value;

Обработка списков (C++2011 variadic templates)

len :: [type] -> int
len x:xs = 1 + len xs
len [] = 0
template <class ...xs>
struct len {};

template <class x, class ...xs>
struct len<x, xs...> { enum { value = 1 + len<xs...>::value }; };

template <>
struct len<> { enum { value = 0 }; };

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

factorial n f = factorial (n - 1) (n * f)
factorial 0 f = f
template <uint64_t n>
using Uint = integral_constant<uint64_t, n>;

using Zero = Uint<0>;
using One = Uint<1>;

template <class N, class M>
struct Sub
  : Uint<N::value - M::value> {};

template <class N, class M>
struct Mul
  : Uint<N::value * M::value> {};

template <class N, class F = One>
struct factorial
  : factorial<Sub<N, One>, Mul<F, N>> {};

template <class F>
struct factorial<Zero, F>
  : F {};

Увы, приведённый код не компилируется. Дело в том, что подстановка шаблонов выполняется немедленно, ленивые вычисления для типов не используются, поэтому компилятор не способен «дойти» до template < class F > struct factorial. Следующий вариант также не будет компилироваться:

template <class N, class F = One>
struct factorial
{
  using base = factorial<Sub<N, One>, Mul<F, N>>;
  using type = typename base::type;
  enum { value = base::value };
};

Впрочем, обойти эту проблему можно, сразу включив условие в определение:

template <class N, class F = One>
struct factorial
  : conditional<N::value != 0,
      factorial<Sub<N, One>, Mul<N, F>>,
      F>::type {};
Вариант 3

Отделение имени от параметров шаблона с помощью вложенных определений. Первый известный мне подобный пример — стандартные аллокаторы C++1998 (rebind). Также применяется в Boost.MPL.

Определение

f x = expr
struct f
{
  template <class T> struct apply { typedef expr type };
};

Применение

f x
typename f::template apply<x>::type

Достоинства:

  • все из варианта 2;
  • удобно передавать в функции высшего порядка (универсальность: достаточно оперировать именами метафункций, являющимися обычными типами, а не шаблонами), в частности, такую метафункцию можно передать саму себе в качестве параметра;

Недостаток: ещё более длинный синтаксис применения (для прямого применения возможно сокращение с помощью template-using, что, впрочем, не спасает при описании выражений внутри метафункций).

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

Обратно к goto

Задача

Дан текст (массив 8-битных символов), получить из него текст с нормализованными переводами строк. Предполагаются возможными следующие варианты оформления переводов строк:

  • 1-символьные: LF (linefeed, символ с кодом 10, ‘\n’, ‘\xa’), CR (carriage return, символ с кодом 13, ‘\r’, ‘\xd’);
  • 2-символьные: CR LF и LF CR.

Любой из этих вариантов должен заменяться на один символ LF.

Данный пример представляется мне не слишком простым и не слишком сложным, более того, можно даже позволить себе считать, что решение вышеприведённой задачи имеет практический смысл :).

Алгоритм

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

  • Все символы, кроме LF и CR просто копируются, в случае же обнаружения CR или LF записываем перевод строки и переходим в состояние «встретился CR» или «встретился LF» соответственно. Если затем встретится не (LF или CR), то записываем встретившийся символ и выходим из этих состояний в исходное состояние.
  • Если встретился тот же символ CR или LF соответственно, то записываем перевод строки и остаёмся в том же состоянии.
  • Если же встретился «перекрёстный» символ LF или CR соответственно, то переходим в новое состояние «комбинация» (отвечающее двухсимвольному переводу строки), выход из которого производится с выводом перевода строки и «другого символа», если таковой встретился, либо возможен возврат в состояния «встретился CR» и «встретился LF» при получении соответствующих символов.
  • Выход происходит, когда исходный массив будет просмотрен целиком.

Впрочем, нетрудно заметить, что по сути состояние «комбинация» совпадает с начальным состоянием «общее», т.е. состояния LF и CR попросту «поглощают» перекрёстный символ. Более подробно и ясно этот алгоритм можно представить в виде схемы конечного автомата.

Конечный автомат

CR LF finite state machine scheme

Схема конечного автомата

Программа

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

Пусть функция, реализующая автомат получает исходный массив в виде двух указателей from и end, первый из которых указывает на начало массива, а последний — на фиктивный элемент за последним, т.е. размер массива равен (end — from). Пусть результат записывается в массив, на начало которого указывает последний параметр функции — указатель to. Пусть функция возвращает указатель на элемент, следующий за последним записываемым (таким образом, будет легко определить, сколько символов было записано).

Итак, получим

char* normalize(const char *from, const char *end, char *to)
{
  char ch;
Begin_:
  if (from == end) return to;
  ch = *from++;
  if (ch == '\n') goto LF_;
  if (ch == '\r') goto CR_;
  *to++ = ch;
  goto Begin_;

LF_:
  *to++ = '\n';
  if (from == end) return to;
  ch = *from++;
  if (ch == '\r') goto Begin_;
  if (ch == '\n') goto LF_;
  *to++ = ch;
  goto Begin_;

CR_:
  *to++ = '\n';
  if (from == end) return to;
  ch = *from++;
  if (ch == '\n') goto Begin_;
  if (ch == '\r') goto CR_;
  *to++ = ch;
  goto Begin_;
}

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

char* normalize(const char *from, const char *end, char *to)
{
  char ch;
  while (from != end)
  {
    switch (ch = *from++)
    {
    case '\n':
      do
      {
        *to++ = '\n';
        if (from == end)
          return to;
        ch = *from++;
      } while (ch == '\n');
      if (ch != '\r')
        *to++ = ch;
      break;

    case '\r':
      do
      {
        *to++ = '\n';
        if (from == end)
          return to;
        ch = *from++;
      } while (ch == '\r');
      if (ch != '\n')
        *to++ = ch;
      break;

  default:
    *to++ = ch;
  }
  return to;
}

Switch-case

Стандартный «структурный» подход к переводу конечных автоматов на язык программирования заключается во введении перечислимого типа, задающего набор возможных состояний автомата, переменной state, хранящей текущее состояние явно, и кодированию действий автомата в зависимости от текущего состояния в виде конструкции switch (state) { case … }, выполняемой в цикле.

Часто разница в производительности результирующего машинного кода варианта на основе switch-case не отличается заметно от варианта на goto (банальное чтение символа из потока ввода может оказаться достаточно дорогой операцией на фоне затрат на выполнение switch-case), однако у варианта switch-case есть ещё то практическое преимущество, что такой автомат позволяет более гибкое использование (здесь для простоты не будем рассматривать сопроцедуры как вариант реализации автомата) — явно хранимое состояние можно сохранить и потом запустить автомат «со старого места».

В данной задаче можно себе представить загрузку символов в буфер (массив), обработка и… повторная загрузка и обработка, но автомат нельзя каждый раз запускать из состояния Begin, т.к. 2-символьный перевод строки может оказаться разбит между двумя загрузками и в результате интерпретирован как два перевода, а не один. Т.е. автомат должен «помнить» своё последнее состояние.

Однако, никто не запрещает нам реализовать гибридный вариант, который я и приведу здесь.

struct NewLineNormalizer
{
  enum State { Begin, LF, CR };
  State state;

  explicit NewLineNormalizer(State st = Begin)
    : state(st) {}

  char* operator()(const char *from, const char *end, char *to)
  {
    if (from == end)
      return to;

    char ch;
    switch (state)
    {
    case LF: goto LF_;
    case CR: goto CR_;
    default:; // Begin falls-through here
    }

    do
    {
      switch (ch = *from++)
      {
      case '\n':
        do
        {
          *to++ = '\n';
          if (from == end)
            return state = LF, to;
        LF_:
          ch = *from++;
        } while (ch == '\n');
        if (ch != '\r')
          *to++ = ch;
        break;

      case '\r':
        do
        {
          *to++ = '\n';
          if (from == end)
            return state = CR, to;
        CR_:
          ch = *from++;
        } while (ch == '\r');
        if (ch != '\n')
          *to++ = ch;
        break;

      default:
        *to++ = ch;
      }
    } while (from != end);

    state = Begin;
    return to;
  }
};

Ниже UPD 19.04.2015.

Switch-case + goto

Хотя выше и приведены варианты, «избавленные» от goto, выглядят они явно запутаннее, чем goto, отвечающие схеме автомата. В последнем случае («гибридном варианте») всё равно использовался goto для «заскока» в нужное место ранее остановленного автомата. Поэтому я решил добавить ещё один «гибридный вариант» — уже на основе оригинального goto-варианта.

class NewLineNormalizer
{
  enum { Begin, LF, CR } state;

public:
  // Конструктор устанавливает начальное состояние.
  NewLineNormalizer() : state(Other) {}

  // Оператор "скобки" позволяет вызывать объект данного класса как функцию.
  // Возвращает указатель на символ, следующий за последним записанным символом.
  char* operator()(const char *from, const char *end, char *to)
  {
    char input;
    switch (state)
    {
    Begin_:
    case Begin:
      if (from == end) // выход
        return state = Begin, to;
      input = *from++;
      if (input == '\n') goto LF_;
      if (input == '\r') goto CR_;
      *to++ = input;
      goto Begin_;

    LF_:
      *to++ = '\n';
    case LF:
      if (from == end) // выход
        return state = LF, to;
      input = *from++;
      if (input == '\r') goto Begin_;
      if (input == '\n') goto LF_;
      *to++ = input;
      goto Begin_;

    CR_:
      *to++ = '\n';
    case CR:
      if (from == end) // выход
        return state = CR, to;
      input = *from++;
      if (input == '\n') goto Begin_;
      if (input == '\r') goto CR_;
      *to++ = input;
      goto Begin_;

    default:
      assert(false);
      return to;
    }
  }
};

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

В этой записи я решил представить алгоритм, придуманный мной под впечатлением от распределённых (distributed) алгоритмов. Алгоритм строит (субоптимальную) правильную вершинную раскраску неориентированного графа. Алгоритм довольно прост и, возможно, был в том или ином виде представлен в литературе.

Алгоритм

  1. Построить первоначальную раскраску (необязательно правильную). Эвристика: построить раскраску в два цвета поиском в ширину (чётный фронт → цвет «0», нечётный фронт → цвет «1»).
  2. Пока есть конфликты, производить разрешение конфликтов.

Под конфликтами понимается наличие одинаково окрашенных вершин-соседей. Проверка «есть конфликты» перебирает все вершины и составляет список C троек (вершина v, количество конфликтов n, минимальный неконфликтный цвет c), отвечающий конфликтующим вершинам (n = количество соседей v, имеющих тот же цвет, что и v). Если вершину v перекрасить в цвет c, то конфликт будет разрешён (при условии сохранения прежних цветов соседями v). Если конфликтов нет, то список C будет пуст, и алгоритм закончит работу — будет получена правильная раскраска вершин графа.

«Разрешение конфликтов» производится следующим образом.

  1. Завести множество соседей перекрашенных вершин A := Ø.
  2. Упорядочить список конфликтов C по убыванию числа конфликтов n.
  3. Для каждой тройки (v, n, c) из C, если v не принадлежит A, то:
    1. Назначить вершине v цвет c.
    2. Пополнить A соседями v.

Свойства алгоритма

По завершении работы алгоритма получаем правильную вершинную раскраску. Необходимо доказать, что работа алгоритма завершится. Строгого доказательства я выполнять не стал. Однако, нетрудно заметить, что всего конфликтов может быть не больше, чем всего вершин в графе (N), за одно разрешение конфликтов разрешается как минимум один конфликт и не добавляется новых, поэтому достаточно N итераций цикла. Одна итерация цикла требует O(N + M) действий (для каждой вершины проверить всех её соседей). Отсюда можно сделать вывод, что алгоритм имеет сложность O(N(N + M)) ≤ O(N3).

Реализация

Была написана реализация алгоритма на C++, включающая набор простых графов для проверки. Для всех простых тестов алгоритм даёт оптимальную раскраску.

Графы-5: паросочетания

Определения

Паросочетанием называется подмножество рёбер графа, содержащее только взаимонеинцидентные рёбра. Любая вершина графа инцидентна не более, чем одному из рёбер паросочетания. Рёбра, входящие в паросочетания, и вершины, инцидентные этим рёбрам, называются насыщенными.

Максимальным называется паросочетание, к которому нельзя добавить ни одного ребра из множества рёбер графа, или, что то же самое, максимальное паросочетание не является подмножеством никакого другого паросочетания.

Наибольшим называется паросочетание, содержащее наибольшее число рёбер. Наибольшее паросочетание является максимальным, обратное неверно.

Полным называется паросочетание, насыщающее все вершины графа. Если полное паросочетание существует, то любое наибольшее паросочетание будет также полным паросочетанием.

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

Алгоритмы

Введённые определения позволяют поставить следующие задачи.

  1. Построить некоторое наибольшее паросочетание.
  2. Проверить, существует ли полное паросочетание.
  3. Построить оптимальное паросочетание для заданных весов рёбер.

Для простоты ограничимся рассмотрением двудольных графов.

Двудольным графом называется граф, вершины которого можно разделить на два непересекающихся множества X и Y таких, что все рёбра графа будут иметь вид (x, y), где x принадлежит X, а yY. Множества X и Y называются долями графа.

Добавим к графу две вспомогательные вершины s и t. Соединим s с каждой вершиной из X, и каждую вершину из Y соединим с вершиной t. Зададим всем рёбрам единичную пропускную способность и будем искать целочисленный максимальный поток (максимальный 0-1-поток) из s в t на построенной таким образом сети. Рёбра с единичным потоком составят некоторое из наибольших паросочетаний. Поэтому задачу 1 можно решить с помощью алгоритма построения максимального потока. Построив же некоторое наибольшее паросочетание не сложно проверить его на полноту, решив и задачу 2.

Кроме алгоритмов построения максимального потока существует адаптация алгоритма Форда-Фалкерсона к задаче построения наибольшего паросочетания, известная как алгоритм Хопкрофта-Карпа. Для решения задачи 2, а именно определения отсутствия полного паросочетания без построения наибольшего паросочетания, также существуют специальные алгоритмы (например, алгоритм Куна). Для решения задачи 3 существует алгоритм Куна-Манкреса, также известный как венгерский алгоритм.

Пример реализации

Возьмём часть кода до строки 506 и поместим её в файл Graph4.h. Для использования алгоритмов построения максимального потока в качестве алгоритмов поиска наибольшего паросочетания нам понадобится три элемента: а) функция, формирующая матрицу пропускных способностей, б) функция-обёртка алгоритма построения максимального потока, в) функция, «извлекающая» паросочетание из матрицы потока.

Первая и третья функции названы prepareBipartiteForMaxFlow и matchingFromFlow соответственно. Размеры долей графа передаются как константы X и Y. Для представления паросочетания используется вектор индексов, отображающий индекс вершины из доли X (от 0 до X — 1) в индекс вершины из доли Y (от 0 до Y), при этом фиктивное значение Y используется, чтобы показать, что вершина не насыщена.

Собственно привязка осуществляется тривиальной функцией fullMatching, которая использует вспомогательную структуру-шаблон FullMatching<class MaxFlowWrapper>, функцию-член find которого и вызывает (так сделано из-за того, что в C++ нельзя определять специализации функций-шаблонов).

Класс, переданный в качестве параметра MaxFlowWrapper, уже «переадресует» вызов конкретному алгоритму построения максимального потока или решает задачу самостоятельно.

Определены следующие классы, подходящие в качестве параметра MaxFlowWrapper: три обёртки алгоритмов из Graph4.h EdmondsKarpMaxFlow, PreflowPushingMaxFlow и PushRelabelMaxFlow, а также реализация алгоритма Хопкрофта-Карпа HopcroftKarp.

Итак, код реализации и теста: codepad.