Category Archives: Lessons

О вычислении матричной экспоненты

1. Экспонента квадратной матрицы \mathbf A формально определяется как сумма ряда:

e^{\mathbf A} \triangleq \mathbf I + \sum\limits_{i=1}^{\infty} \dfrac{\mathbf A^i}{i!},

где \mathbf I есть единичная матрица соответствующего размера.

Из этого определения несложно получить ряд свойств экспоненты, некоторые из которых перечислены ниже:

  1. Экспонента нулевой матрицы есть единичная матрица.
  2. Транспонирование (а также эрмитово сопряжение) перестановочно с экспоненцированием: e^{\mathbf A^\top} = (e^{\mathbf A})^\top.
  3. Если \mathbf A симметрична (эрмитова), то и e^{\mathbf A} симметрична (эрмитова).
  4. Если \mathbf {A B} = \mathbf {B A}, то e^{\mathbf A} \mathbf B = \mathbf B e^{\mathbf A} и e^{\mathbf A} e^{\mathbf B} = e^{\mathbf B} e^{\mathbf A} = e^{\mathbf A + \mathbf B}. В частности, e^{a \mathbf A} e^{b \mathbf A} = e^{(a + b)\mathbf A}.
  5. Матрица e^{\mathbf A} всегда обратима: e^{\mathbf A} e^{-\mathbf A} = \mathbf I.
  6. Для любой обратимой матрицы \mathbf P верно e^{\mathbf{P A P}^{-1}} = \mathbf P e^{\mathbf A} \mathbf P^{-1}.
  7. Экспонента диагональной матрицы \mathbf A = \mathrm{diag}(a_{ii}) также является диагональной матрицей e^{\mathbf A} = \mathrm{diag}(e^{a_{ii}}).
  8. Последнее свойство обобщается на блочно-диагональные матрицы. Пусть \mathbf A = \mathrm{diag}(\mathbf A_k), где каждая \mathbf A_k есть квадратная матрица. Тогда экспонента также является блочно-диагональной матрицей: e^{\mathbf A} = \mathrm{diag}(e^{\mathbf A_k}).

Таким образом, если матрица \mathbf A диагонализируема, т.е. если найдётся такая матрица \mathbf P, что \mathbf A = \mathbf {P D P}^{-1}, где \mathbf D = \mathrm{diag}(\lambda_i) — диагональная матрица (матрица собственных значений), то после выполнения диагонализации (eigendecomposition) \mathbf A вычислить её экспоненту не составит труда.

К сожалению, не все матрицы диагонализируемы, хотя для матриц над полем комплексных чисел множество недиагонализируемых матриц есть множество меры нуль. Поэтому, кстати, даже экспоненту вещественнозначной матрицы может быть удобнее вычислять в комплексных числах. Те, кто хорошо учил линейную алгебру, вспомнят о «жордановой форме», которая существует для любой матрицы и представляет собой как раз блочно-диагональную матрицу, составленную из блоков, называемых «жордановы клетки». Это означает, что если у нас есть представление \mathbf A = \mathbf{P J P}^{-1}, где \mathbf J = \mathrm{diag}(\mathbf J_k), и если мы умеем вычислять экспоненту жордановой клетки, то мы опять же без проблем вычислим экспоненту исходной матрицы. Ниже я коснусь вопроса вычисления экспоненты жордановой клетки.

Особняком стоит свойство, известное как «формула Якоби». Его можно применять, например, для проверки качества вычисления матричной экспоненты:

\det e^{\mathbf A} = e^{\mathrm{tr}\,\mathbf A}.

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

  • Компенсационное суммирование. Например, алгоритм Ноймайера.
  • Вместо вычисления очередной степени \mathbf A^i домножением, я вычислял степень каждый раз «заново» (до 2\log_2{i} + 1 матричных умножений).
  • Единичная матрица добавляется в конце (опционально: не нужно её добавлять, если вычисляется e^{\mathbf A} - \mathbf I ).
  • Суммирование продолжается до тех пор, пока очередное слагаемое не обнулится из-за слишком большого делителя-факториала (на двойной точности обычно достаточно несколько десятков слагаемых).

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

3. Экспонента жордановой клетки. Жорданова клетка \mathbf J (в случае комплексных чисел) есть квадратная матрица, на главной диагонали которой стоит некое число \lambda (любое значение, но одно во всех позициях диагонали), на первой наддиагонали стоят единицы, а все прочие элементы — нули:

\mathbf J = \begin{bmatrix}  \lambda & 1 & 0 & \ldots & 0 \\  0 & \lambda & 1 & \ldots & 0 \\  & & \cdots & & \\  0 & 0 & \ldots & \lambda & 1 \\  0 & 0 & \ldots & 0 & \lambda  \end{bmatrix}.

Заметим, что \mathbf J = \mathbf \Lambda + \mathbf N, где \mathbf \Lambda — шаровая матрица (диагональная с одинаковыми значениями на диагонали), а \mathbf N — нильпотентна, поскольку \mathbf N^n = \mathbf 0 (здесь и ниже n — размер \mathbf J). Более того, легко заметить, что домножение \mathbf N^i на \mathbf N попросту сдвигает единичную наддиагональ на одну позицию вверх. Поэтому имеем

e^{\mathbf N} = \mathbf I + \sum\limits_{i=1}^{n-1}\dfrac{\mathbf N^i}{i!} = \begin{bmatrix}  1 & 1 & 1/2 & 1/6 & \ldots & 1/(n - 1)! \\  0 & 1 & 1 & 1/2 & \ldots & 1/(n-2)! \\  & & & \cdots & & \\  0 & 0 & \ldots & 0 & 1 & 1 \\  0 & 0 & \ldots & 0 & 0 & 1  \end{bmatrix}.

Ну и по свойствам экспоненты и шаровых матриц имеем

e^{\mathbf J} = e^{\mathbf \Lambda + \mathbf N} = e^{\mathbf \Lambda} e^{\mathbf N} = e^\lambda e^{\mathbf N}.

В случае наличия скалярного множителя t:

e^{\mathbf J t} = e^{\mathbf \Lambda t + \mathbf N t} = e^{\lambda t} e^{\mathbf N t},

где

e^{\mathbf N t} = \begin{bmatrix}  1 & t & t^2/2 & t^3/6 & \ldots & t^{n-1}/(n - 1)! \\  0 & 1 & t & t^2/2 & \ldots & t^{n-2}/(n-2)! \\  & & & \cdots & & \\  0 & 0 & \ldots & 0 & 1 & t \\  0 & 0 & \ldots & 0 & 0 & 1  \end{bmatrix}.

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

4. О вычислении матричной экспоненты (и иных функций, определяемых как сумма ряда) с помощью комбинирования разложения Шура, рекуррентных соотношений Сильвестра и суммировании ряда см. Davies, Higham. A Schur-Parlett Algorithm for Computing Matrix Functions.

Реклама

Об аналитическом решении линейных ОДУ

1. Пусть дана однородная система (задача Коши):

\dot{\mathbf x} = \mathbf A \mathbf x, \quad \mathbf x(t_0) = \mathbf x_0,

где \mathbf x(t) \in \mathbb R^n — фазовый вектор системы, \mathbf A — постоянная n\times n-матрица над полем \mathbb R. Цель — найти функцию \mathbf x(t), t \geqslant t_0, удовлетворяющую указанной системе. Запись \dot{\mathbf x} обозначает производную \mathbf x'(t).

Как и в случае n = 1 решение можно получить через экспоненту.
В нашем случае — матричную экспоненту, определяемую формально как сумму ряда. Вычисление матричной экспоненты, особенно для больших n, есть отдельная большая тема, которой в данной статье я не намерен касаться подробно. Ниже приведены формулы для прямого вычисления матричной экспоненты в случае n = 2.

Итак, сразу решение:

\mathbf x(t) = e^{\mathbf A(t - t_0)} \mathbf x_0.

Заметим, что, как и в скалярном случае, мы можем продифференцировать функцию e^{\mathbf A t} по t как матрицу покомпонентно, получив:

(e^{\mathbf A t})' = e^{\mathbf A t} \mathbf A = \mathbf A e^{\mathbf A t}.

\mathbf A, \mathbf A^{-1} и e^{\mathbf A t} являются перестановочными матрицами, что легко показать, раскрыв экспоненту в сумму ряда.

2. Пусть теперь дана неоднородная система:

\dot{\mathbf x} = \mathbf A \mathbf x + \mathbf f, \quad \mathbf x(t_0) = \mathbf x_0,

где \mathbf f \colon \mathbb R \to \mathbb R^n. Выполнив замену

\mathbf z(t) = e^{\mathbf A(t_0 - t)} \mathbf x(t),

исключим фазовую переменную из правой части уравнения:

\dot{\mathbf z}(t) = e^{\mathbf A(t_0 - t)} \mathbf f(t), \quad \mathbf z(t_0) = \mathbf x_0,

откуда легко получить окончательную формулу для решения исходной неоднородной системы:

\mathbf x(t) = e^{\mathbf A(t - t_0)} \left( \mathbf x_0 + \int_{t_0}^t {e^{\mathbf A(t_0 - \tau)} \mathbf f(\tau)} \,d\tau\right).

3. Пусть \mathbf f(t) \equiv \mathbf u = const. Пусть также \det\mathbf A \neq 0. Тогда имеем первообразную:

(e^{\mathbf A t} \mathbf A^{-1})' = e^{\mathbf A t}.

После подстановки первообразной в интеграл формула решения из предыдущего пункта обретает вид:

\mathbf x(t) = e^{\mathbf A(t - t_0)} (\mathbf x_0 + \mathbf A^{-1} \mathbf u) - \mathbf A^{-1} \mathbf u.

4. Предположим теперь, что матрица \mathbf A является функцией времени. Оставим также функцию \mathbf f. Можно воспользоваться тем же приёмом, что и в п.2, введя замену:

\mathbf z(t) = e^{\int_t^{t_0} {\mathbf A(\tau)}\,d\tau} \mathbf x(t).

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

\forall t_1, t_2 \in [t_0, t] \left(\,\mathbf A(t_1) \mathbf A(t_2) = \mathbf A(t_2) \mathbf A(t_1)\,\right).

Тогда можно записать решение системы в форме:

\mathbf x(t) = e^{\int_{t_0}^t {\mathbf A(\tau)}\,d\tau} \left( \mathbf x_0 + \int_{t_0}^t {e^{\int^{t_0}_\tau {\mathbf A(s)}\,ds} \mathbf f(\tau)} d\tau\right).

5. Формулы для матричной экспоненты размерности 2. Пусть дана матрица

\mathbf A = \begin{bmatrix} a & b \\ c & d  \end{bmatrix}

Положим D = 4 bc + (a - d)^2. Имеем теперь три случая.

Случай D > 0.
Положим q = \sqrt{D}, C = \cosh\frac{q}{2}, S = \sinh\frac{q}{2}. Тогда

e^{\mathbf A} = \frac{1}{q} e^\frac{a+d}{2}  \begin{bmatrix} qC + (a - d) S & 2bS \\                             2cS & qC - (a - d)S  \end{bmatrix}

Случай D < 0.
Формула та же, что и в случае D > 0, если положить q = \sqrt{-D}, C = \cos\frac{q}{2}, S = \sin\frac{q}{2}.

Случай D = 0.

e^{\mathbf A} = e^\frac{a+d}{2}  \begin{bmatrix} 1 + \frac{a - d}{2} & b \\                             c & 1 - \frac{a - d}{2}  \end{bmatrix}

Для множества частных случаев можно вывести более простые формулы. Например, для диагональной матрицы b = c = 0:

e^{\mathbf A} = \begin{bmatrix}  e^a & 0 \\ 0 & e^d \end{bmatrix}

Для верхнетреугольной матрицы c = 0:

e^{\mathbf A} = \begin{bmatrix}  e^a & b\frac{e^a - e^d}{a - d}   \\ 0 & e^d \end{bmatrix}

Один алгоритм раскраски графа 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]

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

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