Tag Archives: C++

Об использовании 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».

О задаче разложения числа в сумму

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

Дано: x \in \mathbb R, n \in \mathbb N, \{\underline{x}_i\}_{i=1}^n \subset \mathbb R, \{\overline{x}_i\}_{i=1}^n \subset \mathbb R, \forall i \in \overline{1,n}\, (\underline{x}_i \leqslant \overline{x}_i).

Найти \{x_i\}_{i=1}^n такие, что \sum_{i=1}^n{x_i} = x и \forall i \in \overline{1,n}\, (\underline{x}_i \leqslant x_i \leqslant \overline{x}_i), или определить, что таких x_i не существует.

Задача легко решается, если догадаться сделать замену переменных:

y_i = x_i - \underline{x}_i, \quad \overline{y}_i = \overline{x}_i - \underline{x}_i.

y = \sum\limits_{i=1}^n{y_i} = \sum\limits_{i=1}^n{(x_i - \underline{x}_i)} = x - \sum\limits_{i=1}^n{\underline{x}_i}.

Она позволяет перейти к более простой постановке:

Дано: y \in \mathbb R, n \in \mathbb N, \{\overline{y}_i\}_{i=1}^n \subset \mathbb R, \forall i \in \overline{1,n}\, (0 \leqslant \overline{y}_i).

Найти \{y_i\}_{i=1}^n такие, что \sum_{i=1}^n{y_i} = y и \forall i \in \overline{1,n}\, (0 \leqslant y_i \leqslant \overline{y}_i), или определить, что таких y_i не существует.

Сразу видно, что решение существует тогда и только тогда, когда 0 \leqslant y \leqslant \sum_{i=1}^n{\overline{y}_i} и получить его можно простым линейным жадным алгоритмом, вычитая из y на каждом шаге максимально возможный очередной недостаток y_i = \min\{\,y - \sum_{j=1}^{i-1}{y_j}, \overline{y}_i\,\} (естественно, y_1 = \min\{\,y, \overline{y}_1\}).
Если мы быстро исчерпаем y, то последние y_i могут быть просто нулями.

Приведу код решения данной задачи на C++ (аж с ноября не публиковал здесь кода C++!).

#include <cassert>
#include <iterator>
#include <algorithm>
/**
 * @brief Decompose a number into a sum of non-negative terms.
 * @param value the value to be decomposed
 * @param ubound non-negative upper bounds for decomposition terms
 * @param decomp decomposition terms, of the same size as ubound!
 * @return success (whether a decomposition exists)
 */
template 
  <
    class NT,       // number type
    class NTArrayU, // number type array for ubound
    class NTArrayD  // number type array for decomp
  >
bool scalar_sum_decompose
  (
    NT               value,
    NTArrayU const & ubound,
    NTArrayD       & decomp
  )
{
  if (value < NT{})
    return false;
  
  using std::size;
  auto const n = size(ubound);
  assert(n == size(decomp));
  
  NT s {}; // accumulated sum.
  for (auto i = n-n; i < n; ++i)
  {
    assert(!(ubound[i] < NT{}));
    NT const term = std::min<NT>(value - s, ubound[i]);
    decomp[i] = term;
    value -= term;
  }
  
  return value == NT{};
}

(Код «auto i = n-n» есть маленький костыль ввиду нежелания городить извлечение типа размера.)
Теперь несложно написать и обёртку, решающую первую задачу.

#include <vector>
#include <numeric>
#include <functional>
/**
 * @brief Decompose a number into a sum of bounded terms.
 * @param value the value to be decomposed
 * @param lbound lower bounds for decomposition terms
 * @param ubound upper bounds for decomposition terms
 * @param decomp decomposition terms, it must have the same size as ubound and lbound!
 * @return success (whether a decomposition exists)
 */
template 
  <
    class NT,       // number type
    class NTArrayL, // number type array for lbound
    class NTArrayU, // number type array for ubound
    class NTArrayD  // number type array for decomp
  >
bool scalar_sum_decompose
  (
    NT               value,
    NTArrayL const & lbound,
    NTArrayU const & ubound,
    NTArrayD       & decomp
  )
{
  using std::size;
  auto const n = size(ubound);
  assert(size(lbound) == n && size(decomp) == n);
  
  using std::begin;
  using std::end;
  std::transform(
    begin(ubound), end(ubound), begin(lbound),
    begin(decomp), std::minus<>{});
  NT const lsum = std::accumulate(
    begin(lbound), end(lbound), NT{});
  std::vector<NT> buffer(n);
  if (!scalar_sum_decompose(value - lsum, decomp, buffer))
    return false;
  std::transform(
    begin(buffer), end(buffer), begin(lbound),
    begin(decomp), std::plus<>{});
  return true;
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    mass += cm
    return mass, l

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

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

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

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

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

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

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

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

*out = val;

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

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

#include <type_traits>

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

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

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

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

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

В стандартной библиотеке математических функций C++ (унаследованной от C) имеется функция hypot(x, y), возвращающая евклидову длину вектора (x, y).

Зачем особая функция, если можно просто определить «в лоб»?

inline float hypot_v1(float x, float y)
{
  return std::sqrt(x * x + y * y);
}

(Здесь я использую float и только float по причине, о которой будет сказано ниже.)

Причин может быть несколько:

  1. Гарантировать отсутствие переполнения (из-за возведения в квадрат большого числа или сложения больших чисел) в случае, когда конечный результат (длина) представим.
  2. Гарантировать отсутствие влияния исчезновения (из-за возведения в квадрат слишком маленького числа) на результат. В частности, в случае, когда один из аргументов x, y есть нуль, а другой — число, hypot(x, y) должен возвращать модуль другого аргумента, даже если это субнормальное число.
  3. Обеспечивать погрешность лучше 1 ULP.

Легко показать, что hypot_v1 не удовлетворяет всем этим пунктам. Экспериментально установлено, что погрешность hypot_v1 на основной части диапазона, где не возникает ситуация вредоносного переполнения или исчезновения составляет 1 ULP (проверка для y из набора { 0.f, 1e-40f, 1e-30f, 1e-20f, 1e-15f, 1e-6f, 1e-2f, 1.f, 1e+2f, 1e+6f, 1e+15f, 1e+20f, 1e+30f } и для всех возможных значений x).

В активе hypot_v1 сравнительно большая скорость вычисления (у меня получилось, что она примерно в 10 раз быстрее, чем стандартная hypot на g++ -O3).

Можно ли получить какой-то иной компромисс в плане скорости и качества? В Википедии можно наблюдать следующий вариант (на момент написания). Ниже дан исправленный код:

inline float hypot_v2(float x, float y)
{
  x = std::fabs(x);
  y = std::fabs(y);

  const float
    min_ = x < y? x: y,
    max_ = x < y? y: x;

  if (max_ == 0.f || std::isinf(max_))
    return max_;

  const float q = min_ / max_;
  return max_ * std::sqrt(1.f + q*q);
}

Хорош ли этот код? Он страхует нас от патологического переполнения и исчезновения. Что можно сказать о его скорости? В моём эксперименте получилось, что hypot_v2 примерно вдвое медленнее hypot_v1, что можно считать приемлемым. Что можно сказать о погрешности? Экспериментально получено максимальное значение 2 ULP.

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

inline float hypot_v2a(float x, float y)
{
  x = std::fabs(x);
  y = std::fabs(y);

  if (x == 0.f)
    return y;
  if (y == 0.f)
    return x;

  float
    min_ = x < y? x: y,
    max_ = x < y? y: x;

  int min_e, max_e;
  min_ = std::frexp(min_, &min_e);
  max_ = std::frexp(max_, &max_e);

  min_ *= min_;
  max_ *= max_;

  return std::ldexp(
    std::sqrt(std::ldexp(min_, 2*(min_e - max_e)) + max_),
    max_e);
}

Этот вариант действительно обеспечивает погрешность <= 1 ULP (экспериментально). К сожалению, этот код в моём случае работал втрое медленнее стандартной hypot. (Что может говорить о потенциально низком быстродействии кода предложенного мной для перемножения массивов чисел.)

Теперь следует прояснить вопрос об экспериментальном вычислении погрешности. Так как используются числа float (IEEE-754 binary32, 23 двоичных разряда после запятой), и доступны числа double (IEEE-754 binary64, 52 двоичных разряда после запятой), то нетрудно организовать вычисление ближайшего к истинной величине значения в формате float, попросту выполнив вычисление с точностью double и округлив затем результат до float:

inline float hypot_v3(float x, float y)
{
  const double x_ = x, y_ = y;
  return (float)std::sqrt(x_ * x_ + y_ * y_);
}

Кстати, этот вариант также не страдает от переполнения и исчезновения и вычисляется лишь вдвое медленнее hypot_v1, обеспечивая погрешность <= 0.5 ULP. Он плох тем, что не масштабируется на double (хотя можно задействовать софтовое вычисление с точностью IEEE-754 binary128, заплатив за это резким падением производительности), зато очень удобен для определения качества других вариантов, включая стандартную hypot.

В случае g++ результат hypot везде совпал с результатом hypot_v3 (т.е. точность <= 0.5 ULP стандартной hypot обеспечена). Интересно, что для её вычисления используется древний (1993г.) код от Sun Microsystems, довольно зубодробительного вида.

В принципе, альтернативой hypot и hypot_v1 может стать их гибрид, использующий hypot_v1 там, где она работает хорошо (с погрешностью 1 ULP), а в остальных случаях (предположительно, редких) вызывающий hypot:

inline float hypot_v4(float x, float y)
{
  x = std::fabs(x);
  y = std::fabs(y);

  if (1e-19f <= x && x <= 1e+19f &&
      1e-19f <= y && y <= 1e+19f)
    return std::sqrt(x*x + y*y);

  return std::hypot(x, y);
}

Такой код способен работать со скоростью, в среднем близкой к скорости hypot_v1.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

#include <cmath>
#include <iterator>

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

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

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

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

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

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

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

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

  Number significand() const
  {
    return p;
  }

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

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

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

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

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

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

Conar 6: Parser — окончание

Начало

Предыдущая часть

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

Все нижеприведённые функции являются членами класса Parser.

Функция, проверяющая параметр на удовлетворение первому варианту («ключ»):

bool is_direct_key(const std::string &param) const
{
  return mapper.count(param) != 0;
}

Функция, проверяющая параметр на удовлетворение второму варианту («ключ=значение»), возвращает позицию знака «=»:

size_t is_key_equals(const std::string &param) const
{
  const auto eq_pos = param.find('=');
  if (eq_pos == std::string::npos)
    return 0;

  return is_direct_key(param.substr(0, eq_pos))? eq_pos: 0;
}

Функция, проверяющая параметр на удовлетворение третьему варианту («-ключключ…ключ»):

bool is_key_sequence(const std::string &param) const
{
  const auto psz = param.size();
  if (psz < 3 || param[0] != '-')
    return false;

  char key[3] {'-'};
  for (size_t i = 1; i < psz; ++i)
  {
    key[1] = param[i];
    if (mapper.count(key) == 0)
      return false;
  }

  return true;
}

Наконец, функция, проверяющая параметр на удовлетворение четвёртому варианту («ключзначение»), возвращает позицию первого символа значения:

size_t is_key_value(std::string param) const
{
  while (!param.empty())
  {
    param.pop_back();
    if (is_direct_key(param))
      return param.size();
  }

  return 0;
}

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

Используя определённые выше функции можно набросать скелет цикла обработки последовательности параметров:

  Unparsed_parameters unparsed;
  Parameters params(from, to);
  for (auto key = begin(params), params_end = end(params),
            prev_key = params_end;
    key != params_end;)
  {
    // Check all possible variants one by one.
    if (is_direct_key(*key))
    {
      // The only situation in which we may have a sequence
      // of argument parameters after the key.

      prev_key = key;
    }
    else if (const auto eq_pos = is_key_equals(*key))
    {    }
    else if (is_key_sequence(*key))
    {    }
    else if (const auto val_pos = is_key_value(*key))
    {    }
    else
    {
      // Nope, this is not a valid parameter of any form.
      unparsed.emplace_back(move(*key));
      ++key; // next...
    }
  }
  return unparsed;

Если сработало условие is_direct_key(*key), то надо найти следующий ключ, а всё что между передать диапазоном обработчику опции, соответствующей *key. Первая мысль: запоминать последнюю такую позицию (итератор prev_key) и вызывать обработчик, когда был найден правый конец. Недостаток: это надо делать под каждым if’ом, хотелось бы избежать такого дублирования кода…

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

    size_t pos;
    if
    (
      int situation =
        is_direct_key(*key)?          1
      : (pos = is_key_equals(*key))?  2
      : is_key_sequence(*key)?        3
      : (pos = is_key_value(*key))?   4
      : 0 // --> else
    )
    {   }
    else
    {
      // Nope, this is not a valid parameter of any form.
      unparsed.emplace_back(move(*key));
      ++key; // next...
    }

Теперь я могу избежать дублирования кода, а необходимая информация о случае хранится в situation и pos. Хотя нет, есть ещё случай, когда ключ — последний, дальше до конца идут только его аргументы… Продублировать вызов обработчика после цикла? Думаю, лучше ситуацию «конец параметров» тоже занести в situation:

  for (auto key = begin(params), params_end = end(params),
            prev_key = params_end;;)
  {
    size_t pos;
    if
    (
      int situation =
        key == params_end?            1
      : is_direct_key(*key)?          2
      : (pos = is_key_equals(*key))?  3
      : is_key_sequence(*key)?        4
      : (pos = is_key_value(*key))?   5
      : 0 // --> else
    )

Интересно, как можно было бы записать то же самое чище и без дублирования кода…

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

template <class Key>
Handler& handler_of(const Key &key)
{
  return std::get<0>(options[mapper.at(key)]);
}

Теперь я могу, наконец, записать, что делается под if’ом:

  if (prev_key != params_end)
  {
    for (auto last = handler_of(*prev_key)
                     (prev_key + 1, key);
         last != key; ++last)
    { unparsed.emplace_back(move(*last)); }
    prev_key = params_end;
  }

  switch (situation)
  {
  case 1:
    return unparsed;

  case 2:
    prev_key = key;
    break;

  case 3: case 5:
    {
      const int has_equals = situation == 3;
      auto k = key->substr(0, pos);
      key->erase(0, pos + has_equals);
      if (handler_of(k)(key, key + 1) == key)
        unparsed.emplace_back(move(*key));
    }
    break;

  case 4:
    for (size_t i = 1, sz = key->size(); i < sz; ++i)
    {
      char k[3] { '-', (*key)[i] };
      handler_of(k)(params_end, params_end);
    }
  }

Красотой этот код не отличается, увы. И надо написать для него тесты.

В результате тестирования была исправлена ветка else:

    else
    {
      // Nope, this is not a valid parameter of any form.
      if (prev_key == params_end)
        unparsed.emplace_back(move(*key));
    }

При тестировании со строковыми параметрами обнаружилась старая ошибка: оператор >> читает «словами», однако в случае строковых параметров их надо не «читать» (с пробельными символами-разделителями, да), а просто «отдавать» целиком. Для этого я заменю явное чтение с помощью >> на вызов новой функции impl::read_value:

/// Tries to read a generic val from string str.
template <class S, class T>
inline bool read_value
  (std::istringstream &reader, S &&str, T &val)
{
  reader.str(std::forward<S>(str));
  return reader >> val;
}

/// Forwards a string instead of reading it.
inline bool read_value
  (std::istringstream&, std::string &&str, std::string &val)
{
  val = std::move(str);
  return true;
}

Забавно, что после этой замены сработал тест value_out (код ошибки 7) — из-за того, что теперь «чтение» забирает строки с помощью move, а не копирует их.

Вместе с тестами таки вылезло за пределы 1000 строк. В итоге, я не решился использовать C++17 (только static_assert без сообщения, а можно было декомпозиционное определение и string_view, например).

Код данного варианта целиком.

О затенении перегрузок функций

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

Рассмотрим простенький код: определим оператор вывода вектора в поток.

#include <iostream>
#include <iterator>
#include <vector>
using namespace std;

ostream& operator<<(ostream &os, const vector<int> &v)
{
	for (auto x: v)
		os << x << ' ';
	return os;
}

int main()
{
	cout << vector<int>{ 1, 2, 3, 4 } << endl;
	return 0;
}

Он работает так, как и предполагается, никаких проблем нет. Range-for — удобная штука: до его появления обход контейнера с помощью итераторов выглядел страшненько. Вместо явного цикла можно было воспользоваться возможностями STL: алгоритмом copy и адаптером потока ostream_iterator.

#include <iostream>
#include <iterator>
#include <vector>
#include <algorithm>
using namespace std;

ostream& operator<<(ostream &os, const vector<int> &v)
{
	ostream_iterator<int> oi(os, " ");
	copy(begin(v), end(v), oi);
	return os;
}

int main()
{
	cout << vector<int>{ 1, 2, 3, 4 } << endl;
	return 0;
}

Этот вариант работает точно так же, как и предыдущий. Всё в порядке.

Попробуем обобщить его, например, на вектор векторов…

#include <iostream>
#include <iterator>
#include <vector>
#include <algorithm>
using namespace std;

ostream& operator<<(ostream &os, const vector<int> &v)
{
	ostream_iterator<int> oi(os, " ");
	copy(begin(v), end(v), oi);
	return os;
}

ostream& operator<<
	(ostream &os, const vector<vector<int>> &m)
{
	ostream_iterator<vector<int>> oi(os, "\n");
	copy(begin(m), end(m), oi);
	return os;
}

int main()
{
	using VI = vector<int>;
	vector<VI> m
	{
		VI{ 1, 2, 3, 4 },
		VI{ 5, 6, 7, 8 },
		VI{ 9, 0, 1, 2 },
	};
	
	cout << m << endl;
	return 0;
}

И вот она, ошибка… Этот код уже не компилируется. Как часто бывает в таких случаях, компилятор может высыпать ворох маловразумительных ошибок. Проблема в операторе <<. Точнее, в том факте, что в том же пространстве имён std, где определены ostream_iterator и vector, определён и набор стандартных перегрузок operator<<, который «затеняет» operator<< из глобального пространства имён.

Простейший способ избавиться от ошибки — поместить свою версию operator<< также в пространство имён std. Например, такой код уже благополучно компилируется и работает правильно:

#include <iostream>
#include <iterator>
#include <vector>
#include <algorithm>
using namespace std;

namespace std
{
	ostream& operator<<(ostream &os, const vector<int> &v)
	{
		ostream_iterator<int> oi(os, " ");
		copy(begin(v), end(v), oi);
		return os;
	}
}

ostream& operator<<
	(ostream &os, const vector<vector<int>> &m)
{
	ostream_iterator<vector<int>> oi(os, "\n");
	copy(begin(m), end(m), oi);
	return os;
}

int main()
{
	using VI = vector<int>;
	vector<VI> m
	{
		VI{ 1, 2, 3, 4 },
		VI{ 5, 6, 7, 8 },
		VI{ 9, 0, 1, 2 },
	};
	
	cout << m << endl;
	return 0;
}

Я намеренно не стал помещать второй operator<< внутрь namespace std, чтобы показать, что здесь проблема была с первым operator<<, который должен вызываться некоей функцией, определённой в std.

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

Следующий пример определяет рекурсивную версию operator<<, которую тоже приходится поместить в namespace std:

#include <iostream>
#include <iterator>
#include <vector>
#include <algorithm>
using namespace std;

namespace std
{
	template <class T>
	ostream& operator<<(ostream &os, const vector<T> &m)
	{
		ostream_iterator<T> oi(os.put('{'), " ");
		copy(begin(m), end(m), oi);
		return os.put('}');
	}
}

int main()
{
	using VI = vector<int>;
	vector<VI> m
	{
		VI{ 1, 2, 3, 4 },
		VI{ 5, 6, 7, 8 },
		VI{ 9, 0, 1, 2 },
	};
	
	cout << m << endl;
	return 0;
}

Вообще, «по умолчанию» считается плохим тоном вносить собственные определения в пространство имён std (есть исключения, например, частные специализации hash).

Для операций с собственными типами, определёнными в собственном пространстве имён, возможно положиться на ADL, поместив перегрузки операторов в то же самое пространство имён. В вышеприведённых примерах использовался std::vector, а определение операций для стандартных классов приводит к опасности коллизии разных определений для одних и тех же шаблонных классов, появившихся в разных местах проекта. Лично мне бы хотелось иметь в C++ «strong typedef», тем более, что для него даже новых ключевых слов вводить не надо:

namespace my
{
	using Data = new vector<int>;
	// Теперь Data != vector<int> и
	// определён в пространстве имён my.
	ostream& operator<<(ostream &os, const Data &d);
}

Мечты-мечты.

Conar 5: структура данных Parser

Начало

Предыдущая часть

Займёмся классом Parser.

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

Обработчик будем приводить к типу Handler, объявленному следующим образом:

using Parameters = std::vector<std::string>;
using PIt = Parameters::iterator;
using Handler = std::function<PIt(PIt, PIt)>;

Теперь определим собственно структуры данных, содержащиеся в объекте Parser. Mapper отображает ключи в индексы вектора Options.

using Option = std::tuple<Handler, Option_info>;
using Options = std::vector<Option>;
using Mapper = std::map<std::string, std::size_t>;
Options options;
Mapper mapper;

Напишем код, добавляющий опцию (тройку «обработчик, ключи, инфо«).

/// Add an option description.
/// Use functions flag, value and seq to make an option.
template <class Opt>
Parser& operator()(Opt option)
{
  using namespace std;
  static_assert(tuple_size<Opt>::value == 3);

  // Add the option to options.
  const auto index = options.size();
  options.emplace_back(move(get<0>(option)),
                       move(get<2>(option)));

  // Register keys of the option.
  for (auto &key: get<1>(option))
    mapper.emplace(key_prep(move(key)), index);

  return *this;
}

Функция key_prep добавляет знаки «-» в начало строки:

static std::string key_prep(std::string &&key)
{
  if (!key.empty() && std::isalnum(key[0]))
  {
    if (key.size() == 1)
      key.insert(0, 1, '-');
    else
      key.insert(0, "--");
  }

  return move(key);
}

Теперь можно заполнить тело функции help, которую можно использовать и в целях тестирования. Подумав, я решил убрать параметр «ширина строки», поскольку это не очень нужное усложнение. Реальную ширину строки консоли нельзя получить стандартными методами, а насильная вставка перевода строки может приводить к появлению избыточных пустых строк. Функция получилась длинная, но простая:

/// Construct a help message.
std::string help(
  std::size_t key_column_width = 20,
  const char *line_sep = "\n\n") const
{
  using namespace std;
  ostringstream write;
  write.setf(ios::left);

  // Option index -> first key (iterator).
  vector<Mapper::const_iterator> key_by_index(options.size());
  for (auto p = mapper.end(), b = mapper.begin();;)
  {
    --p;
    key_by_index.at(p->second) = p;
    if (p == b)
      break;
  }

  // Make the text.
  for (auto p = mapper.begin(), e = mapper.end(); p != e; ++p)
  {
    if (p->first.size() < key_column_width)
    {
      write.width(key_column_width);
      write << p->first;
    }
    else
    {
      write.width(0);
      write << p->first << '\n';
    }

    const auto pivot = key_by_index[p->second];
    if (pivot == p)
      write << get<1>(options[p->second]);
    else
      write << ("See " + pivot->first + ".");
    write << line_sep;
  }

  return write.str();
}

Самую сложную часть парсера — собственно разбор списка параметров — из-за недостатка времени в данный момент я оставил для следующей части.

Код данного варианта целиком.

Продолжение следует…