Category Archives: Programming

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Пример:

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

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

Конкатенация и транспонирование tuple

В недавно вышедшей Boost 1.66 beta 1 присутствует новая библиотека метапрограммирования Mp11. Я начал читать отсюда и обратил внимание на два примера.

Первый — конкатенация произвольного набора кортежей типов. Мне показалось очень сложным представленное решение Эрика Ниблера, поскольку «на глаз» было видно, что стандартными средствами это сделать довольно легко, ЕСЛИ ограничиться поддержкой std::tuple. Естественно, я попробовал написать свою версию, которую и привожу ниже.

Вначале конкатенация двух tuple:

template <class A, class B>
struct Tcat;

template <class... T, class... S>
struct Tcat<std::tuple<T...>, std::tuple<S...>>
{
  using type = std::tuple<T..., S...>;
};

template <class A, class B>
using tcat = typename Tcat<A, B>::type;

Теперь не составит труда определить конкатенацию произвольного количества tuple рекурсивно:

template <class... A> struct Tcats;

template <class... T>
struct Tcats<std::tuple<T...>>
{
  using type = std::tuple<T...>;
};

template <class A, class B, class... Other>
struct Tcats<A, B, Other...>
{
  using type = typename Tcats<tcat<A, B>, Other...>::type;
};

template <class... T>
using tcats = typename Tcats<T...>::type;

Второй — transform бинарной операцией двух кортежей (zip, если в терминах функциональных языков программирования). Конкретный случай для двух кортежей не столь интересен. Легко видеть, что при обобщении на произвольное число кортежей мы получаем ни что иное как транспозицию матрицы (F — применяемая операция):


<T11, T12, T13, ..., T1m>,
<T21, T22, T23, ..., T2m>,
...
<Tn1, Tn2, Tn3, ..., Tnm>
--->
F<T11, T21, ..., Tn1>,
F<T12, T22, ..., Tn2>,
F<T13, T23, ..., Tn3>,
...
F<T1m, T2m, ..., Tnm>

Насколько просто сделать такую штуку на шаблонах C++? Результат, который у меня получился, мне не особенно понравился. Думаю, можно это сделать лучше, однако, не факт, что это «лучше» будет означать «короче».

Итак, метафункция, порождающая из столбца исходной матрицы строку новой матрицы с префиксом F (n-арной операцией):

template <class I,
  template <class...> class F, class... T>
struct Tcol
{
  using type = F<std::tuple_element_t<I::value, T>...>;
};

template <class I,
  template <class...> class F, class... T>
using tcol = typename Tcol<I, F, T...>::type;

Здесь I — тип, определяющий зависимое имя value, раскрываемое в целочисленную константу — индекс столбца.

Теперь надо применить tcol к диапазону индексов столбцов. Так как готовых iota, map и bind у меня нет, я написал рекурсивную реализацию аналога «for» «в лоб»:

template <class S, class T,
  template <class...> class F, class... Other>
struct Tfor;

template <class Int, Int s,
  template <class...> class F, class... Other>
struct Tfor<
  std::integral_constant<Int, s>,
  std::integral_constant<Int, s>,
  F, Other...>
{
  using type = std::tuple<>;
};

template <class Int, Int s, Int t,
  template <class...> class F, class... Other>
struct Tfor<
  std::integral_constant<Int, s>,
  std::integral_constant<Int, t>,
  F, Other...>
{
  static_assert(s < t, "");
  using S  = std::integral_constant<Int, s>;
  using S1 = std::integral_constant<Int, s+1>;
  using T  = std::integral_constant<Int, t>;
  
  using type = tcat<std::tuple<F<S, Other...>>,
    typename Tfor<S1, T, F, Other...>::type>;
};

template <class S, class T,
  template <class...> class F, class... Other>
using tfor = typename Tfor<S, T, F, Other...>::type;

Здесь S и T — типы, задающие начало и конец диапазона (::value — соответствующие целочисленные константы), пак Other — исходная матрица типов.

Теперь записать обобщённый zip не составит труда:

template <
  template <class...> class F,
  class... T>
struct Tzip;

template <
  template <class...> class F>
struct Tzip<F> { using type = std::tuple<>; };

template <
  template <class...> class F,
  class T0, class... T>
struct Tzip<F, T0, T...>
{
  template <class I, class... Y>
  using tcolF = tcol<I, F, Y...>;
  using type = tfor<
    std::integral_constant<std::size_t, 
      0>,
    std::integral_constant<std::size_t,
      std::tuple_size<T0>::value>,
    tcolF, T0, T...>;
};

template <
  template <class...> class F,
  class... T>
using tzip = typename Tzip<F, T...>::type;

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

Тестировочный код:

template <class... T>
struct Tmax;

template <class T>
struct Tmax<T> { using type = T; };

template <class T0, class... T>
struct Tmax<T0, T...>
{
  using tail_max = typename Tmax<T...>::type;
  using type = std::conditional_t<
    (T0::value < tail_max::value), tail_max, T0>;
};

template <class... T>
using tmax = typename Tmax<T...>::type;

template <class T>
using Sizeof =
  std::integral_constant<std::size_t, sizeof(T)>;

template <class... T>
using tmax_sizeof = tmax<Sizeof<T>...>;

int main()
{
  using namespace std;
  // Concat.
  {
    using t0 = tuple<int, int*>;
    using t1 = tuple<float>;
    using t2 = tuple<char*, char**, void*, void**>;
    using t3 = tuple<>;
    using ex = tuple<
      int, int*, float, char*, char**, void*, void**>;
    static_assert(is_same<tcats<t0, t1, t2, t3>, ex>::value);
  }
  
  // Generic zip.
  {
    using r0 = tuple<char, short, int, long>;
    using r1 = tuple<unsigned, void*, void**, void***>;
    using exr = tuple<
      tuple<char, unsigned>,
      tuple<short, void*>,
      tuple<int, void**>,
      tuple<long, void***>>;
    //print_t(tzip<tuple, r0, r1>{});
    static_assert(is_same_v<exr, tzip<tuple, r0, r1>>);
    
    using t0 = tuple<int[11], int[12], int[13]>;
    using t1 = tuple<int[21], int[22], int[23]>;
    using t2 = tuple<int[31], int[32], int[33]>;
    using t3 = tuple<int[41], int[42], int[43]>;
    using ext = tuple<
      tuple<int[11], int[21], int[31], int[41]>,
      tuple<int[12], int[22], int[32], int[42]>,
      tuple<int[13], int[23], int[33], int[43]>>;
    static_assert(
      is_same_v<ext, tzip<tuple, t0, t1, t2, t3>>);
    
    using exsz = tuple<
      Sizeof<int[41]>, Sizeof<int[42]>, Sizeof<int[43]>>;
    static_assert(
      is_same_v<exsz, tzip<tmax_sizeof, t0, t1, t2, t3>>);
  }
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    mass += cm
    return mass, l

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

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

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

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

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

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

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

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

*out = val;

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

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

#include <type_traits>

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

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

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

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

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

Пытаясь сделать вариант hypot(x, y), вычисляющий правильно округлённый результат для всех пар чисел x и y, я обнаружил, что используемая для проверки hypot_v3 может быть вовсе не так хороша, как я полагал.

В первой заметке я написал, что она обеспечивает погрешность <= 0.5 ULP, что неточно, так как происходит два округления: при вычислении корня в формате double и затем при конвертировании double во float. Технически погрешность <= 0.5 * (1 + DBL_EPSILON) ULP, эта маленькая добавка (DBL_EPSILON == 2.22045e-16) теоретически может сказаться в редких случаях.

Приведу некоторые очевидные свойства hypot(x, y):

1. Коммутативность: hypot(x, y) == hypot(y, x).
2. Чётность: hypot(-x, y) == hypot(x, -y) == hypot(x, y).
3. Округление: если |x| + |y| == max(|x|, |y|), то hypot(x, y) == max(|x|, |y|).

Данные свойства позволяют оценить количество "разных" пар x, y, дающих нетривиальный результат.

Всего имеется неотрицательных чисел float: 255·223 = 2139095040 (включая +бесконечность). Из них субнормальных 223-1 = 8388607 (с нулём), нормальных — 2130706433. Для субнормальных имеем 8388607·4194303 = 35184359505921 комбинаций. Для нормальных — 2130706433·224 = 35747322059030528 комбинаций. Всего 35782506418536449 ≈ 3.58e+16 комбинаций. Если считать, что DBL_EPSILON является хорошей оценкой вероятности ошибки, то умножив её на число комбинаций получим округлённо 8 — матожидание числа неверно округлённых пар. Считанные единицы, на которые можно надеяться натолкнуться только при переборе хотя бы одной восьмой всех возможных сочетаний значений x и y…

Увы, на практике всё несколько хуже. Для нормальных чисел. И много-много хуже для субнормальных! Действительно, оценка погрешности предполагает нормальную форму числа, при которой ULP имеет фиксированный (в относительной погрешности) вес (2-23). В субнормальных числах последняя цифра имеет больший вес вплоть до 100% (когда это единственная единица — в наименьшем представимом положительном числе). Исходя из этого наблюдения, можно попробовать добавить в hypot_v3 коррекцию порядка:

inline float hypot_v3a(float x, float y)
{
  const double x_ = x, y_ = y,
    s = x_ * x_ + y_ * y_;

  if (s >= 0x1p-251)
    return (float)std::sqrt(s);

  return 0x1p-126f * (float)std::sqrt(0x1p+252 * s);
}

Эта небольшая поправка вроде бы вообще не должна влиять на цифры множителя, однако на практике резко уменьшает число дефектных результатов в субнормальной области (на момент написания я не наткнулся ни на одну подозрительную пару из субнормальных чисел). Среди нормальных подозрение вызывают, например, пары (0.0003162f, 1.661635309e-7f), (0.01f, 0.0001590774482f) и (1e+15f, 4.605317338e+15f).

Итак, я заменил проверочную hypot_v3 на hypot_v3a. Мы знаем, что результат sqrt(x*x + y*y) после коррекции порядка отличается от правильно округлённой гипотенузы не более, чем на вес последней цифры. Здесь надо сделать замечание по поводу того, что именно я под этим подразумеваю. Представление чисел в формате IEEE-754 обладает удобным свойством: их побитовым представлением можно оперировать как целым числом. В случае float 22 младших бита есть двоичные цифры множителя после запятой. Перед запятой подразумевается 1, если число нормальное, и 0, если субнормальное. За 22 битами цифр следует 8 бит порядка. Значение порядка 0 соответствует нулю и субнормальным числам. Значения 1—254 соответствуют нормальным числам (множители от 2-126 до 2127). Порядок 255 соответствует бесконечности (нулевые младшие 22 бита) и нечислам (ненулевые).

В моём случае получается, что модуль разности результата hypot_3a и результата вычисления корня одинарной точности, интерпретируемых как целые числа в соответствии с их двоичным представлением, не превосходит единицы. А значит, получив неточный результат z, я могу рассмотреть два соседних числа, добавив или вычтя 1 из двоичного представления z как из целого числа (то же, что std::nextafter(z, INFINITY), std::nextafter(z, 0.f), но стандартная реализация nextafter тоже не блещет скоростью — она предполагает обработку различных случаев, которых у меня не может быть). Итак, если удастся оценить погрешность результата, то помимо вычисленного значения должно быть достаточно проверить следующее (добавить единицу как к целому) и предыдущее (вычесть единицу как из целого) числа — какое-то из них может оказаться правильным ответом.

В начале, впрочем, я приведу улучшенную версию hypot_v2b, которая избегает работы с субнормальными числами.

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

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

  if (max_ <= 0x1p-126f)
  {
    invmul = 0x1p-126f;  // 0x1p-23 <= max_ <= 1
    mul = 0x1p+126f;     // 0x1p-23 <= min_ *
  }
  else if (max_ <= 0x1p-63f)
  {
    invmul = 0x1p-87f;   // 0x1p-39 <= max_ <= 0x1p+24
    mul = 0x1p+87f;      // 0x1p-63 <= min_ *
  }
  else if (max_ <= 1.f)
  {
    invmul = 0x1p-24f;   // 0x1p-39 <= max_ <= 0x1p+24
    mul = 0x1p+24f;      // 0x1p-63 <= min_ *
  }
  else if (max_ <= 0x1p+63f)
  {
    invmul = 0x1p+39f;   // 0x1p-39 <= max_ <= 0x1p+24
    mul = 0x1p-39f;      // 0x1p-63 <= min_ *
  }
  else
  {
    invmul = 0x1p+102f;  // 0x1p-38 <= max_ <  0x1p+26
    mul = 0x1p-102f;     // 0x1p-63 <= min_ *
  }

  max_ *= mul;
  min_ *= mul;

  if (max_ + min_ == max_)
    return invmul * max_;
  // * here.
  return invmul * std::sqrt(max_* max_ + min_* min_);
}

Пусть есть два соседних числа в плавающей точке a и b = a(1-\varepsilon). Пусть истинное значение корня находится между ними в точке z(t) = (1-t)a + t b. Пограничное значение соответствует t = 1/2. Обозначим разности известных квадратов: \delta_1(t) = z(t)^2 - a^2, \delta_2(t) = z(t)^2 - b^2. Дискриминантом назовём величину

d(t) = 2(\delta_1(t) + \delta_2(t)) + (a - b)^2.

Заметим, что дискриминант обращается в нуль при t = 1/2, меньше нуля при t > 1/2 и больше нуля при t < 1/2. Итак, b является правильно округлённым z(t), если d(t) < 0, либо d(t) = 0 и двоичное представление b оканчивается на нуль (округление «серединок» к чётным).

Аналогично, выбрав b = a(1 + \varepsilon), получим, что d(t) больше нуля при t > 1/2.

Теперь воспользуемся double для того, чтобы написать «reference»-версию с оценкой погрешности на основе hypot_v2d.

namespace float_bits
{
  using ifloat = std::int32_t;
  using ufloat = std::uint32_t;
  static_assert(sizeof(ufloat) == sizeof(float));

  inline ufloat to_bits(const float x)
  {
    return reinterpret_cast<const ufloat&>(x);
  }

  inline float to_float(const ufloat bits)
  {
    return reinterpret_cast<const float&>(bits);
  }

  inline bool of_same_sign(float a, float b)
  {
    return ((to_bits(a) ^ to_bits(b)) >> 31) == 0;
  }

  inline float bits_prev(float x)
  {
    return to_float(to_bits(x) - 1);
  }

  inline float bits_next(float x)
  {
    return to_float(to_bits(x) + 1);
  }

  inline ufloat ulp_distance(float a, float b)
  {
    const auto x = to_bits(a), y = to_bits(b);
    return x < y? y - x: x - y;
  }
}

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

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

  if (max_ <= 0x1p-126f)
  {
    invmul = 0x1p-126f;  // 0x1p-23 <= max_ <= 1
    mul = 0x1p+126f;     // 0x1p-23 <= min_
  }
  else if (max_ <= 0x1p-63f)
  {
    invmul = 0x1p-87f;   // 0x1p-39 <= max_ <= 0x1p+24
    mul = 0x1p+87f;      // 0x1p-63 <= min_
  }
  else if (max_ <= 1.f)
  {
    invmul = 0x1p-24f;   // 0x1p-39 <= max_ <= 0x1p+24
    mul = 0x1p+24f;      // 0x1p-63 <= min_
  }
  else if (max_ <= 0x1p+63f)
  {
    invmul = 0x1p+39f;   // 0x1p-39 <= max_ <= 0x1p+24
    mul = 0x1p-39f;      // 0x1p-63 <= min_
  }
  else
  {
    invmul = 0x1p+102f;  // 0x1p-38 <= max_ <  0x1p+26
    mul = 0x1p-102f;     // 0x1p-63 <= min_
  }

  max_ *= mul;
  min_ *= mul;

  if (max_ + min_ == max_)
    return invmul * max_;

  const double a = max_, b = min_, s = a*a + b*b;
  const float
    s_ = (float)s,
    s0 = std::sqrt(s_),
    s1 = float_bits::bits_prev(s0),
    s2 = float_bits::bits_next(s0);

  const double d0 = s0, d1 = s1, d2 = s2,
    q0 = s - d0 * d0,
    q1 = s - d1 * d1,
    q2 = s - d2 * d2,
    aq1 = std::fabs(q1),
    aq2 = std::fabs(q2);

  float res = s0;
  if (aq1 < aq2)
  {
    const auto q1_0 = d1 - d0,
      d = 2. * (q1 + q0) + q1_0 * q1_0;
    if (d < 0. || 
       (d == 0. && (float_bits::to_bits(s1) & 1) == 0))
      res = s1;
  }
  else if (aq2 < aq1)
  {
    const auto q2_0 = d2 - d0,
      d = 2. * (q2 + q0) + q2_0 * q2_0;
    if (d > 0. ||
       (d == 0. && (float_bits::to_bits(s2) & 1) == 0))
      res = s2;
  }

  return invmul * res;
}

Результат hypot_v2e совпадает с результатом hypot_v3a почти везде: проверялись все x из набора { 3.16227766e-4f, 0.f, 1e-40f, 0x1p-126f — 0x1p-127f, 0x1p-126f, 1e-30f, 1e-20f, 1e-15f, 1e-6f, 1e-2f, 1.f, 1e+2f, 1e+6f, 1e+15f, 1e+20f, 1e+30f } для всех неотрицательных представимых чисел y — было найдено всего четыре пары (x, y), для которых результаты не совпали — три из них были приведены выше как «подозрительные».

Я написал аналог hypot_v2e, который оперирует только значениями float. Точное возведение в квадрат выполняется с помощью разрезания множителя числа на две 12-битные половинки (hi и lo). Затем вычисляются слагаемые hi*hi, lo*lo и 2*hi*lo. Результат — суммарно 50 отличающихся от hypot_v3a результатов (для того же перебора x-y), т.е. в среднем 1 ошибка на несколько сот миллионов пар (x, y). Тем не менее, не скажу, что такой результат мне особенно понравился. Хотелось бы почти полного совпадения с hypot_v2e или большей простоты кода…

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

В посте Эксперименты с hypot я допустил существенную неточность: указал ссылку на реализацию std::hypot (double) без проверки реализации std::hypotf (float), в то время как для сравнения использовалась как раз версия одинарной точности. Итак, hypotf.

Оказывается, это, по сути, та же hypot_v3! Я было думал, что hypotf использует ту же методику, что и hypot. В этой мысли меня поддерживало и её низкое быстродействие (в 5 раз медленнее hypot_v3). Видимо, дополнительные проверки оказались столь дороги. Таким образом, с одной стороны, очевидно, что результаты hypot_v3 и hypotf (для чисел) совпадут. С другой стороны, разработчики руководствовались теми же соображениями, что и я. Функция hypot_v3 даёт корректно округлённый результат в одинарной точности, и её можно использовать для экспериментальной проверки прочих вариантов.

Вариант в двойной точности не мог быть создан с помощью того же приёма, так как четверная точность не поддерживается аппаратно на большинстве процессоров, а программная реализация очень медленная и требует подключения специальной библиотеки (хотя она и есть в комплекте gcc). И вот тут возник вопрос: float vs double даёт возможность проверить качество того или иного алгоритма, а если мы его применяем только в double, то мы, по сути, ограничены теоретическим анализом.

Я, наконец, прочёл комментарий в e_hypot.c, который, вроде бы (с кодом не сравнивал), описывает применённый там метод. А метод основан на соображении следующего вида:

«If (assume round-to-nearest) z=x*x+y*y has error less than sqrt(2)/2 ulp, than sqrt(z) has error less than 1 ulp (exercise).»

Т.е. если вычислить правильно округлённую сумму квадратов, то и корень будет правильно округлён (докажите, ага). Конечно, чисто формально можно разложить в ряд:

\sqrt{z+\varepsilon}-\sqrt{z} = \dfrac{\varepsilon}{2\sqrt{z}} + O(\varepsilon^2).

Если привести абсолютную погрешность к относительной, то получим

\dfrac{\sqrt{z(1+\varepsilon)} - \sqrt{z}}{\sqrt{z}} = \sqrt{1+\varepsilon}-1 = \dfrac{\varepsilon}{2} + O(\varepsilon^2).

В первом приближении извлечение корня уполовинивает относительную погрешность! С такой точки зрения, у нас должен получаться правильно округлённый результат даже в том случае, если последняя цифра z неверна. Впрочем, следует учесть погрешность вычисления самого корня. Если считать, что эпсилон = 0.5 ULP, то у нас появляется дополнительная относительная погрешность \varepsilon:

\left(1 + \dfrac{\varepsilon}{2} + O(\varepsilon^2)\right)(1 + \varepsilon) = 1 + \dfrac{3}{2}\varepsilon + O(\varepsilon^2).

Т.е. всё же ошибка в 0.75 ULP. Я выше написал «правильно округлён», но это, увы, не то же самое, что «ошибка меньше 1 ULP». Ошибку меньше 1 ULP мы получили, но для «правильно округлён» требуется ошибка не более 0.5 ULP, чего мы не достигаем. Эксперимент показывает, что, действительно, вычисление вида

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

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

  if (max_ + min_ == max_)
    return max_;

  // max_ != 0
  // min_ >= 0x1p-24 * max_

  if (max_ <= 0x1p-126f)
  {
    invmul = 0x1p-126f;  // 0x1p-23 <= max_ <= 1
    mul = 0x1p+126f;     // 0x1p-23 <= min_
  }
  else if (max_ <= 0x1p-63f)
  {
    invmul = 0x1p-63f;   // 0x1p-63 <= max_ <= 1
    mul = 0x1p+63f;      // 0x1p-87 <= min_
  }
  else if (max_ <= 1.f)
  {
    invmul = 1.f;        // 0x1p-63 <= max_ <= 1
    mul = 1.f;           // 0x1p-87 <= min_
  }
  else if (max_ <= 0x1p+63f)
  {
    invmul = 0x1p+63f;   // 0x1p-63 <= max_ <= 1
    mul = 0x1p-63f;      // 0x1p-87 <= min_
  }
  else
  {
    invmul = 0x1p+127f;  // 0x1p-64 <= max_ <  2
    mul = 0x1p-127f;     // 0x1p-88 <= min_
  }

  max_ *= mul;
  min_ *= mul;

  const double a = max_, b = min_, s = a*a + b*b;
  return invmul * std::sqrt((float)s);
}

даёт ошибку в последней двоичной цифре, хотя и немного реже, чем вариант hypot_v2b. Следовательно, то же поведение можно ожидать от std::hypot: ошибка в последней цифре. Т.е. она немногим лучше прямого вычисления в духе hypot_v2b. В связи с этим возникает вопрос: можно ли более-менее эффективно вычислить действительно правильно округлённую гипотенузу, не прибегая к арифметике удвоенной разрядности. И ещё: как быть с длиной n-мерного вектора? Метод hypot_v2b можно применить к n-мерному вектору: найти максимум по абсолютной величине, выбрать порядок, вычислить сумму промасштабированных квадратов с компенсацией и потом уже взять корень, и можно надеяться на погрешность в одну-две последних цифры.

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

Небольшое дополнение к предыдущему посту. Ввиду низкого быстродействия frexp/ldexp может быть разумно заменить их на выбор множителя (степени двойки) банальными if-else.

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

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

  if (max_ <= 0x1p-126f)
  {
    invmul = 0x1p-126f;
    mul = 0x1p+126f;
  }
  else if (max_ <= 0x1p-63f)
  {
    invmul = 0x1p-63f;
    mul = 0x1p+63f;
  }
  else if (max_ <= 1.f)
  {
    invmul = 1.f;
    mul = 1.f;
  }
  else if (max_ <= 0x1p+63f)
  {
    invmul = 0x1p+63f;
    mul = 0x1p-63f;
  }
  else
  {
    invmul = 0x1p+126f;
    mul = 0x1p-126f;
  }

  max_ *= mul;
  min_ *= mul;

  return invmul * std::sqrt(max_* max_ + min_* min_);
}

Данный код предотвращает некорректное переполнение и исчезновение и обеспечивает погрешность 1 ULP. Конечно, гребёнка if-ов сказывается на скорости работы — у меня получилось, что hypot_v2b примерно втрое медленнее hypot_v1 (но всё же в десять раз быстрее hypot_v2a). Эта оценка довольно условна, так как скорость исполнения ветвлений зависит от множества факторов.

Аналогичный способ можно применить для замены frexp/ldexp при вычислении произведений. Например, если произведение получилось меньше по абсолютной величине наименьшего положительного нормального числа, то пересчитать его, предварительно домножив на 0x1p+63f; если произведение стало бесконечным (переполнение) или приблизилось по величине к наибольшему представимому числу, то пересчитать его, домножив на 0x1p-63f.