Category Archives: Lessons

Об использовании Visual Studio Code для учебных примеров

Visual Studio Code — редактор на основе платформы Electron, который можно использовать как минималистичную IDE. Мне требовалась поддержка C++17 и переносимость всего «на флешке». Что меня приятно удивило, так это «шустрость» VSCode.

Итак, порядок действий, предпринятых мной был таков:

  1. Скачать VSCode в архиве (в моём случае это была версия 1.27.1 для Windows/x64). Распаковать его, допустим, в папку vscode. Внимание! Путь к vscode не должен содержать не-ASCII символы и пробелы.
  2. Создать пустую папку vscode/data. Там будут плагины и общие настройки. При обновлении VSCode можно просто распаковать новую версию и переместить туда свою папку data из старой.
  3. Распаковать сборку MinGW рядом с vscode, допустим, в папку mingw (mingw/bin содержит g++.exe). В качестве таковой я брал сборку отсюда, там уже есть Boost и ряд других библиотек. Я к ним ещё добавляю Eigen и кое-что своё. Добавляю, просто складывая соответствующие папки в mingw/include.
  4. Создать рядом с vscode пустую папку для примеров. Допустим, cpp.
  5. Запустить VSCode (vscode/code.exe) и установить плагин C/C++ for Visual Studio Code от Microsoft (с Microsoft Intellisense).
  6. Открыть в VSCode папку cpp.
  7. Настройки создаются в виде JSON-файлов в папке cpp/.vscode. Ниже я опишу возможное их содержимое. Их можно создать непосредственно в VSCode, заполнить в нём же, и затем перезапустить VSCode.
  8. После того, как файлы настроек созданы, и VSCode перезапущен, следует протестировать его. Создать cpp/test.cpp с произвольным содержимым на C++ (например, helloworld). Terminal/Run build task… или Ctrl+Shift+B должно запускать сборку и завершаться без ошибок. F5 запускает сборку и результат в отладчике (gdb). Заголовочные файлы должны успешно открываться (правой кнопкой мыши и Go to definition). Intellisense должен работать адекватно.
  9. Полученный «пакет» в составе vscode, mingw, cpp можно упаковать в архив и развернуть на другой системе или просто скопировать на внешний носитель и запускать с него.

Итак, файлы настроек.

Плагин C/C++ использует файл c_cpp_properties.json.
Его содержимое в моём случае:

{
    "configurations": [
        {
            "name": "MinGW",
            "intelliSenseMode": "gcc-x64",
            "compilerPath": "${workspaceFolder}\\..\\mingw\\bin\\g++.exe",
            "cppStandard": "c++17",
            "cStandard": "c11",
            "includePath": [
                "${workspaceFolder}",
                "${workspaceFolder}\\..\\mingw\\include"
            ],
            "defines": [
                "__cplusplus=201703L"
            ]
        }
    ],
    "version": 4
}

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

Для сборки требуется создать «задачу сборки». В простейшем случае это компиляция одной единицы трансляции (cpp-файла) сразу в исполняемый файл (в моём случае всегда один и тот же с названием lab.exe).
Описание задач сборки располагается в файле tasks.json. Его содержимое в моём случае:

{
    // See https://go.microsoft.com/fwlink/?LinkId=733558
    // for the documentation about the tasks.json format
    "version": "2.0.0",
    "tasks": [
        {
            "label": "build lab",
            "type": "shell",
            "command": "${workspaceFolder}\\..\\mingw\\bin\\g++",
            "args": [
                "-g", "-Og", // "debug" build
                //"-O3",
                "${file}", 
                "-march=native", "-Wall", "-pedantic",
                "-std=c++17",
                "-olab.exe"
            ],
            "group": {
                "kind": "build",
                "isDefault": true
            }
        }
    ]
}

Наконец, описание запуска отладчика по F5 располагается в launch.json.
Его содержимое в моём случае:

{
    // Use IntelliSense to learn about possible attributes.
    // Hover to view descriptions of existing attributes.
    // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
    "version": "0.2.0",
    "configurations": [
        {
            "name": "(gdb) Launch",
            "type": "cppdbg",
            "request": "launch",
            "program": "${workspaceFolder}\\lab.exe",
            "args": [],
            "stopAtEntry": false,
            "cwd": "${workspaceFolder}",
            "environment": [],
            "externalConsole": true,
            "MIMode": "gdb",
            "miDebuggerPath": "${workspaceFolder}\\..\\mingw\\bin\\gdb.exe",
            "setupCommands": [
                {
                    "description": "Enable pretty-printing for gdb",
                    "text": "-enable-pretty-printing",
                    "ignoreFailures": true
                },
                {
                    "description": "Enable break on exception for gdb",
                    "text": "catch throw",
                    "ignoreFailures": true              
                }
            ],
            "preLaunchTask": "build lab"
        }
    ]
}

Вторая команда catch throw активирует перехват исключений отладчиком gdb.
Вообще, к сожалению, я не могу сказать, что с отладкой всё хорошо: есть ряд странностей. Например, точка останова, поставленная в одном месте, может сработать ниже (как будто поставленная в другом месте — может быть, это можно исправить, если заменить -Og на -O0). Точку останова нельзя поставить внутри шаблона (не срабатывает вообще). Не даёт ставить новые точки останова во время отладки. Если во время отладки нажать паузу, то продолжение может не срабатывать (останется только остановить исполнение).

Для меня проблемы с отладчиком не являются критичными, так как я привык к отладке через проверку условий и протоколирование. Но тем, кто привык к отладчику в «нормальной» Visual Studio такое, конечно, будет не по нраву.

UPD. Если всё это запускать на Windows 10, то в качестве терминала может быть по умолчанию выбран PowerShell. Почему-то он «не понимает» относительные пути через .. и не может запустить компилятор. В этом случае рекомендуется заменить его на cmd.exe (в VSCode можно выбрать терминал).

UPD2. Замена флага оптимизации -Og на -O0 всё же помогает отладчику, по крайней мере, точки останова не «убегают» с тех мест, куда они были поставлены.

Реклама

Граница Парето в 2D: STL-версия

Почти пять лет назад я описал алгоритм вычисления границы Парето множества точек на плоскости.

С тех пор я использовал данную задачу как задачу по программированию на первом курсе (на алгоритмы STL). Всё дело в том, что этот алгоритм очень просто записывается с помощью STL. Требуется применить два алгоритма (sort и unique) и адаптер компаратора (not_fn или самописный, например, переставляющий аргументы местами). Бывают простенькие задачи, к которым невозможно адекватно применить STL без написания своих компонент, а бывает и наоборот. Итак, код.

#include <algorithm>
#include <functional>

template 
  <
    class RanIt, 
    class Comp1, 
    class Comp2
  >
RanIt remove_dominated
  (
    RanIt from,
    RanIt to,
    Comp1 compare1,
    Comp2 compare2
  )
{
  std::sort(from, to, std::not_fn(compare1));
  return std::unique(from, to, std::not_fn(compare2));
}

////////////////////////////
#include <string>
#include <vector>
#include <iostream>
using namespace std;

bool test_remove_dominated()
{
  vector<string> input
  {
    "look",
    "after",
    "the",
    "raven",
    "star"
  };
  
  input.erase(
    remove_dominated(
      input.begin(), input.end(),
      std::less<>{},
      [](auto & a, auto & b)
      {
        return a.size() < b.size();
      }),
    input.end()
    );
    
  sort(input.begin(), input.end());
  return input == vector<string>
    {
      "raven",
      "star",
      "the"
    };
}

int main()
{
  cout << test_remove_dominated() << '\n';
  return 0;
}

В примере два критерия сравнения строк: лексикографический и по длине. Так «raven» доминирует «look» (и лексикографически, и по длине), но не «star».

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

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! Я провёл несколько тестов данного кода с многочленами разного размера и он дал правильный результат, однако я не даю никаких гарантий относительно его корректности и, тем более, реального быстродействия (это я не проверял).