Парсер Хабра: 2021-09-26

суббота, 2 октября 2021 г.

Дизайн и эволюция constexpr в C++

constexpr - одно из самых магических ключевых слов в современном C++. Оно дает возможность создать код, который будет выполнен еще до окончания процесса компиляции, что является абсолютным пределом для быстродействия программ.

У constexpr с каждым годом становится больше возможностей. Сейчас использовать в compile-time вычислениях можно почти всю стандартную библиотеку. Пример вычисления числа до 1000 с наибольшим количеством делителей: ссылка на код.

История constexpr насчитывает долгую историю эволюции с ранних версий C++. Исследуя предложения в стандарт и исходники компиляторов, можно понять, как слой за слоем создавалась эта часть языка, почему именно так она выглядит, как на практике вычисляются constexpr-выражения, какие возможности ждут нас в будущем, а какие - могли бы быть, но не были приняты в Стандарт.

Эта статья подходит как тем, кто еще не знает, что такое constexpr, так и тем, кто уже долгое время его использует.

C++98 и C++03: Сословия среди const-переменных

В C++ в некоторых местах нужно использовать целочисленные константы, значения которых должны быть известны в compile-time. Стандарт разрешает записывать константы в виде несложных выражений, как в этом коде:

enum EPlants {
    APRICOT = 1 << 0,
    LIME = 1 << 1,
    PAPAYA = 1 << 2,
    TOMATO = 1 << 3,
    PEPPER = 1 << 4,
    FRUIT = APRICOT | LIME | PAPAYA,
    VEGETABLE = TOMATO | PEPPER,
};

template<int V> int foo();
int foo6 = foo<1+2+3>();
int foo110 = foo<(1 < 2) ? 10*11 : VEGETABLE>();

int v;
switch (v) {
case 1 + 4 + 7:
case 1 << (5 | sizeof(int)):
case (12 & 15) + PEPPER:
    break;
}

Эти выражения описаны в разделе [expr.const] и называются constant-expression. Они могут содержать только:

Литералы (туда входят целые числа, это интегральные типы)
Значения enum-ов
Параметр шаблона интегрального или enum-типа (напр., значение V из template<int V>)
sizeof-выражение
const-переменные, инициализированные constant-expression - интересный пункт

Все пункты, кроме последнего, неудивительны - они точно известны во время компиляции. Переменные - более интересный случай.

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

В стандартах C++98/03 есть два вида static initialization:

zero-initialization, память заполняется нулями, затем в ходе программы изменяется
initialization with a constant expression, память (если нужна) сразу содержит вычисленное значение

[Note: Все остальные инициализации называются dynamic initialization, мы их не рассматриваем. — end note]

[Note: zero-initialized переменная во время работы программы в свою очередь может быть проинициализирована еще раз "нормальным" значением, это уже будет dynamic initialization (пусть даже до старта main). — end note]

Рассмотрим пример с обоими видами переменных:

int foo() {
    return 13;
}

const int test1 = 1 + 2 + 3 + 4; // initialization with a const. expr.
const int test2 = 15 * test1 + 8; // initialization with a const. expr.
const int test3 = foo() + 5; // zero-initialization
const int test4 = (1 < 2) ? 10 * test3 : 12345; // zero-initialization
const int test5 = (1 > 2) ? 10 * test3 : 12345; // initialization with a const. expr.

Переменные test1, test2, test5 можно использовать как параметр шаблона, значение справа от case в switch и т.д., а test3 и test4 - нельзя.

Как видно из требований к constant-expression и из примера, есть транзитивность - если какая-то составная часть выражения не является constant-expression, то и само выражение не является constant-expression. При этом рассматриваются только те части выражения, которые реально вычисляются - поэтому test4 и test5 попадают в разные группы.

Если нигде не берется адрес constant-expression переменной, то скомпилированной программе разрешено не резервировать для нее память, поэтому "заставим" ее это сделать. Выведем значения переменных и их адреса:

int main() {
    std::cout << test1 << std::endl;
    std::cout << test2 << std::endl;
    std::cout << test3 << std::endl;
    std::cout << test4 << std::endl;
    std::cout << test5 << std::endl;

    std::cout << &test1 << std::endl;
    std::cout << &test2 << std::endl;
    std::cout << &test3 << std::endl;
    std::cout << &test4 << std::endl;
    std::cout << &test5 << std::endl;
}

izaron@izaron:~/cpp$ clang++ --std=c++98 a.cpp 
izaron@izaron:~/cpp$ ./a.out 
10
158
18
180
12345
0x402004
0x402008
0x404198
0x40419c
0x40200c

Скомпилируем объектный файл и посмотрим на таблицу символов:

izaron@izaron:~/cpp$ clang++ --std=c++98 a.cpp -c
izaron@izaron:~/cpp$ objdump -t -C a.o

a.o:     file format elf64-x86-64

SYMBOL TABLE:
0000000000000000 l    df *ABS*  0000000000000000 a.cpp
0000000000000080 l     F .text.startup  0000000000000015 _GLOBAL__sub_I_a.cpp
0000000000000000 l     O .rodata        0000000000000004 test1
0000000000000004 l     O .rodata        0000000000000004 test2
0000000000000004 l     O .bss   0000000000000004 test3
0000000000000008 l     O .bss   0000000000000004 test4
0000000000000008 l     O .rodata        0000000000000004 test5

В конкретной версии компилятора под конкретную архитектуру, в конкретной программе zero-initialized переменные попали в секцию .bss, остальные в секцию .rodata.

Загрузчик операционной системы перед запуском загружает программу так, что секция .rodata оказывается в сегменте с read-only режимом, который защищен от записи на уровне ОС.

Попробуем через const_cast изменить данные по адресу переменных. Стандарт дает уклончивый ответ на тему того, когда запись в результат const_cast-а может спровоцировать Undefined Behaviour. По крайней мере, этого не происходит, если мы убираем const с объекта/указателя на объект, который не является фундаментально константным изначально. Т.е. надо различать физическую константность и логическую.

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

    const int& ref = testX;
    const_cast<int&>(ref) = 13; // OK for test3, test4; SEGV for test1, test2, test3
    std::cout << ref << std::endl;

Таким образом, одни const-переменные "константнее" других. Насколько известно, в то время не было простого способа проверить или как-то проконтролировать, что переменная была initialized with a const. expr.

0-∞: Вычислитель констант в компиляторе

Чтобы представлять, как константые выражения вычисляются во время компиляции, нужно понимать устройство компилятора.

Компиляторы идейно похожи друг на друга, я опишу процесс вычисления на основе компилятора Clang/LLVM. Я скопировал базовую информацию про этот компилятор из своей прошлой статьи:

Clang и LLVM

Про само устройство Clang и LLVM написано уже много статей. На хабре я бы посоветовал эту статью, чтобы понять их краткую историю и общую схему.

Количество стадий компиляций зависит от того, кто объясняет устройство компилятора. Анатомия компилятора многоуровнева и на самом абстрактном уровне выглядит так, что есть три разные программы:

Front-end: переводит исходник из C/C++/Ada/Rust/Haskell/... в LLVM IR - особое промежуточное представление. Фронтендом для C-like языков является Clang.
Middle-end: LLVM IR оптимизируется в зависимости от настроек.
Back-end: LLVM IR переводится в машинный код под нужную платформу - x86/Arm/PowerPC/...

Для простых языков реально написать компилятор под 1000 строк и получить всю мощь фреймворка LLVM - для этого нужно реализовать фронтенд. Также можно использовать lex/yacc - готовые синтаксические парсеры.

На менее абстрактном уровне находится фронтенд Clang, который выполняет такие действия (не рассматривая препроцессор и прочие "микро"-шаги):

Лексический анализ: перевод символов в токены, например []() { return 13 + 37; } преобразуются в (l_square) (r_square) (l_paren) (r_paren) (l_brace) (return) (numeric_constant:13) (plus) (numeric_constant:37) (semi) (r_brace).
Синтаксический анализ: создание AST (Abstract Syntax Tree), то есть перевод токенов из предыдущего пункта в вид (lambda-expr (body (return-expr (plus-expr (number 13) (number 37))))).
Кодогенерация: создание LLVM IR по данному AST.

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

"Микросервис", который занимается вычислением константных выражений (от самых простых в C++98 до сложных в C++23), назовем условно вычислитель констант.

Если согласно Стандарту в каком-то месте кода ожидается константное выражение; и выражение, которое там находится, действительно удовлетворяет требованиям на константное выражение, то Clang должен "не отходя от кассы" в 100% случаев уметь вычислять его.

Ограничения на константные выражения на протяжении лет постоянно смягчались, а вычислитель констант Clang соответсвенно постоянно усложнялся, вплоть до управления моделью памяти.

Есть старая документация 9-летней давности, описывающая вычисление констант для C++98/03. Так как константные выражения тогда были очень простыми, они выполнялись с помощью обычного constant folding через анализ синтаксического дерева (AST). Так как в синтаксическом дереве все арифметические выражения уже разобраны в виде под-деревьев, то вычисление константы - просто элементарный обход под-дерева.

Исходник вычислителя констант находится в lib/AST/ExprConstant.cpp и на момент написания статьи разросся до почти 16 тысяч строк. С годами он научился интерпретировать много всего, например циклы (EvaluateLoopBody), и всё это на синтаксическом дереве.

У константных выражений есть важное отличие от кода, выполняющегося в рантайме - они обязаны не допускать undefined behaviour. Если вычислитель констант наткнется на UB, компиляция провалится.

c.cpp:15:19: error: constexpr variable 'foo' must be initialized by a constant expression
    constexpr int foo = 13 + 2147483647;
                  ^     ~~~~~~~~~~~~~~~

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

c.cpp:15:18: warning: overflow in expression; result is -2147483636 with type 'int' [-Winteger-overflow]
    int foo = 13 + 2147483647;
                 ^

2003: Макросы не нужны

Изменения в стандарт происходят через предложения.

Где находятся предложения и из чего они состоят?

Все предложения в Стандарт находятся на open-std.org. Практически все написаны понятно - чаще всего есть:

Краткий обзор области со ссылками на разделы Стандарта;
Текущие проблемы;
Предлагаемое решение проблем;
Предлагаемые изменения в текст Стандарта;
Ссылки на предложения-предшественники и прошлые ревизии предложения;
В особо продвинутых предложениях - ссылки на реализацию в форке компилятора. В тех, что я видел, реализацию делали в форке clang-a.

По ссылкам на предыдущие предложения можно отследить эволюцию для каждого куска C++.

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

Принять участие в эволюции C++ может любой - для русскоязычных экспертов есть stdcpp.ru.

Предложение от 2003 года [N1521] Generalized Constant Expressions указывает на проблему того, что если часть выражения вычисляется с использованием вызова метода, то выражение не является constant-expression. Это заставляет злоупотреблять макросами, если нужно получить более-менее сложное константное выражение:

#define SQUARE(X) ((X) * (X))
inline int square(int x) { return x * x; }
// ^^^ определение макроса и метода
square(9)
std::numeric_limits<int>::max()
// ^^^ невозможны в составе constant-expression
SQUARE(9)
INT_MAX
// ^^^ теоретически могут быть в составе constant-expression

Поэтому предлагается ввести понятие constant-valued методов, которых будет разрешено использовать в constant-expression. Метод считается constant-valued, если это inline-метод, нерекурсивный, возвращающий не void, и его тело состоит из единственного выражения вида return expr;, где после подстановки аргументов (куда также идут constant-expression) получился бы constant-expression.

[Note: Забегая вперед, термин constant-valued не прижился end note].

int square(int x) { return x * x; } // constant-valued
long long_max(int x) { return 2147483647; } // constant-valued
int abs(int x) { return x < 0 ? -x : x; } // constant-valued
int next(int x) { return ++x; } // NOT constant-valued

Таким образом, все переменные test1-5 из прошлого раздела становились бы "фундаментально" константными, без изменения в коде.

Предложение считает, что можно пойти еще дальше и надо рассмотреть вариант, что такой код тоже должен компилироваться:

struct cayley {
    const int value;
    cayley(int a, int b)
        : value(square(a) + square(b)) {}
    operator int() const { return value; }
};

std::bitset<cayley(98, -23)> s; // eq. to bitset<10133>

Потому что переменная value "фундаментально константная" - инициализировалась в конструкторе через constant-expression с двумя вызовами constant-valued метода. То есть, согласно общей логике предложения, для этот код можно свести примерно к такому (вынос переменных и методов вне структуры):

// имитация вызова конструктора cayley::cayley(98, -23) и operator int()
const int cayley_98_m23_value = square(98) + square(-23);

int cayley_98_m23_operator_int() {
    return cayley_98_m23_value;
}

// создание битсета
std::bitset<cayley_98_m23_operator_int()> s; // eq. to bitset<10133>

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

[Note: Однако предложения не существуют в отрыве от компиляторов - если предложение нереально реализовать в разумный срок, его вряд ли примут end note].

Как с переменными, программист не может проконтролировать, что метод является constant-valued.

2006-2007: Тайное становится явным

К счастью, через 3 года в следующих ревизиях этого предложения ([N2235]) признали, что слишком много неявности это плохо. В список проблем добавили невозможность контроля за инициализацией:

struct S {
    static const int size;
};

const int limit = 2 * S::size; // dynamic initialization
const int S::size = 256; // constant expression initialization
const int z = std::numeric_limits<int>::max(); // dynamic initialization

По задумке программиста, limit должен был быть инициализирован constant-expression, но этого не происходит, так как S::size определён "слишком поздно", после limit. Если бы была возможность запросить нужный тип инициализации, компилятор выдал бы ошибку.

Аналогично с методами. constant-valued методы переименовали в constant-expression методы. Требования к ним остались те же, но теперь, чтобы их было возможно использовать в constant-expression, их необходимо объявлять с ключевым словом constexpr. Зато компиляция будет падать, если тело метода не является правильным return expr;.

Также компиляция упадет, если constexpr-метод принципиально не сможет быть использован в constant-expression при любых аргументах, с ошибкой constexpr function never produces a constant expression. Это надо для того, чтобы программист был точно уверен в том, что метод является потенциально используемым в constant-expression.

Некоторые методы стандартной библиотеки (напр. из std::numeric_limits) предлагается пометить constexpr, так как они удовлетворяют требованиям на него.

Переменные или члены класса также можно объявлять constexpr, тогда компиляция упадет, если переменная инициализируется не через constant-expression.

На тот момент решили сохранить совместимость нового слова с переменными, неявно инициализированными через constant-expression без ключевого слова constexpr, то есть этот код работал (забегая вперед, сейчас этот код с --std=c++11 не компилируется, возможно этот код так никогда и не начал работать):

const double mass = 9.8;
constexpr double energy = mass * square(56.6); // OK, хотя mass не объявлен с constexpr
extern const int side;
constexpr int area = square(side); // error: square(side) is not a constant expression

constant-expression-конструкторы для пользовательских классов также легализовали. Этот конструктор должен иметь пустое тело и инициализировать члены constexpr-expression-ами, если пользователь создает constexpr-объект данного класса.

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

Пример класса с constexpr-ами из предложения:

struct complex {
    constexpr complex(double r, double i) : re(r), im(i) { }

    constexpr double real() { return re; }
    constexpr double imag() { return im; }

private:
    double re;
    double im;
};

constexpr complex I(0, 1); // OK -- literal complex

Такие объекты как I в предложении назвали user-defined-literals. "Литерал" - это нечто вроде базовой сущности в C++. Так же, как "простые" литералы (числа, символы и т.д.) сразу подставляются в ассемблерные команды, а строковые литералы лежат в секции подобной .rodata, так и пользовательские литералы занимают где-нибудь там своё место.

Теперь constexpr-переменными могли являться не только числа и перечисления, а литеральные типы, которых ввели в этом предложении (пока еще без reference type). Литеральный тип - такой, который может быть передан в constexpr-метод, и/или изменен и/или возвращен из нее. Это достаточно простые типы, чтобы компиляторы могли его поддерживать в вычислителе констант.

Ключевое слово constexpr стало спецификатором, нужным компилятору - примерно как override в классах. После обсуждения предложения, для этого ключевого слова решили не делать новый storage class (хотя было бы законно), новый type qualifier, и также его решили не разрешать использовать для аргументов методов, чтобы не переусложнять правила перегрузки методов.

2007: Первые constexpr для структур данных

В этом году вышло предложение [N2349] Constant Expressions in the Standard Library, где пометили как constexpr некоторые методы и константы, а также некоторые методы контейнеров, например:

template<size_t N>
class bitset {
    // ...
    constexpr bitset();
    constexpr bitset(unsigned long);
    // ...
    constexpr size_t size();
    // ...
    constexpr bool operator[](size_t) const;
};

Конструкторы инициализируют члены класса через constant-expression, остальные методы внутри себя имеют return expr;, подходящий под текущие ограничения.

Больше половины предложений про constexpr за все годы - про то, чтобы пометить как constexpr какие-нибудь методы из стандартной библиотеки. Они возникают сразу после очередного витка эволюции constexpr и почти всегда не очень интересны - будем рассматривать изменения в сам язык.

2008: Рекурсивные constexpr-методы

constexpr-методы изначально не предполагалось разрешать делать рекурсивными по причине отсутствия убедительных доводов в наличие рекурсии, но затем это ограничение убрали, что отметили в [N2826] Issues with Constexpr.

constexpr unsigned int factorial( unsigned int n )
{ return n==0 ? 1 : n * factorial( n-1 ); }

У компилятора существует некий предел вложенности вызовов (в clang это 512 вложенных вызовов), при превышении он откажется считать выражение.

Подобные пределы также есть, например, для инстанцирования шаблонов (если мы бы считали в compile-time через шаблоны, а не через constexpr-методы)

2010: "const T&" как аргументы в constexpr-методах

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

template< class T >
constexpr const T& max( const T& a, const T& b ); // не скомпилируется

constexpr pair(); // можно поставить constexpr
pair(const T1& x, const T2& y); // нельзя поставить constexpr

Предложение [N3039] Constexpr functions with const reference parameters (a summary) разрешает константные ссылки в аргументах и как возвращаемое значение.

Это опасное изменение: до этого вычислитель констант имел дело с простыми выражениями и constexpr-переменными (объект литерального класса - по сути набор из constexpr-переменных); но введение ссылок пробивает "четвертую стену", потому что это понятие относится к модели памяти, которой в вычислителе нет.

В общем случае работа со ссылками и указателями в constant-expression превращает C++-компилятор в C++-интерпретатор, поэтому накладываются различные ограничения.

Если вычислитель констант может обработать метод с аргументом типа T, то типа const T& тоже сможет - если будет "представлять себе", что для этого аргумента создается "временный объект".

Более-менее сложная работа со ссылками и попытки что-либо поломать не скомпилируются

template<typename T> constexpr T self(const T& a) { return *(&a); }
template<typename T> constexpr const T* self_ptr(const T& a) { return &a; }
template<typename T> constexpr const T& self_ref(const T& a) { return *(&a); }
template<typename T> constexpr const T& near_ref(const T& a) { return *(&a + 1); }

constexpr auto test1 = self(123); // OK
constexpr auto test2 = self_ptr(123); // FAIL, pointer to temporary is not a constant expression
constexpr auto test3 = self_ref(123); // OK
constexpr auto tets4 = near_ref(123); // FAIL, read of dereferenced one-past-the-end pointer is not allowed in a constant expression

2011: static_assert в constexpr-методах

Предложение [N3268] static_assert and list-initialization in constexpr functions вводит возможность писать "статические" объявления, не вляющие на результат работы метода: typedef, using, static_assert . Это небольшое развинчивание гаек для constexpr-методов.

2012: (Почти) любой код в constexpr-функциях

В 2012 году произошел большой рывок вперед с предложением [N3444] Relaxing syntactic constraints on constexpr functions. Есть множество простых методов, которых желательно уметь вычислять в compile-time, например степень a^n :

// Compute a to the power of n
int pow(int a, int n) {
  if (n < 0)
    throw std::range_error("negative exponent for integer power");
  if (n == 0)
    return 1;
  int sqrt = pow(a, n/2);
  int result = sqrt * sqrt;
  if (n % 2)
    return result * a;
  return result;
}

Однако программистам фокусничать писать в функциональном стиле (убирать локальные переменные и if-ы), чтобы сделать ее constexpr-вариант:

constexpr int pow_helper(int a, int n, int sqrt) {
  return sqrt * sqrt * ((n % 2) ? a : 1);
}
// Compute a to the power of n
constexpr int pow(int a, int n) {
  return (n < 0) ? throw std::range_error("negative exponent for integer power") :
         (n == 0) ? 1 : pow_helper(a, n, pow(a, n/2));
}

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

Так как в константных выражениях нельзя изменять переменные, циклы (for/while/do/range-based for) заведомо невозможно использовать
switch и goto запрещены, чтобы вычислитель констант не моделировал сложные потоки управления
Как при старых ограничениях, для метода теоретически должен существовать набор аргументов, при которых ее можно использовать в константных выражениях. Иначе считается, что метод помечен constexpr ошибочно, и компиляция упадет с constexpr function never produces a constant expression.

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

[Note: таких переменных не будет сильно много из-за жесткого ограничения на "глубину вызовов" end note].

В методах можно объявлять статические переменные. Они могут иметь не-литеральный тип (чтобы, например, возвращать ссылки на них из метода; сами ссылки - тип как раз литеральный), но у них не должно быть dynamic initialization (т.е. должна быть хотя бы zero-initialization) и нетривиального деструктора. Предложение приводит пример, где это могло бы быть полезно (получение ссылки на нужный объект в compile-time):

constexpr mutex &get_mutex(bool which) {
  static mutex m1, m2; // non-const, non-literal, ok
  if (which)
    return m1;
  else
    return m2;
}

Также разрешили объявлять типы (class, enum и т.д.) и возвращать void.

2013: (Почти) любой код в constexpr-функциях ver 2.0 Mutable Edition

Однако Комитет решил, что поддержка циклов (хотя бы for) в constexpr-методах это must have. В 2013 году вышла поправленная версия предложения [N3597] Relaxing constraints on constexpr functions.

Для реализации "constexpr for"-а рассматривалось четыре варианта.

Самым далеким от "остального C++" варианта было создание совершенно новой конструкции для итераций, который подходил бы для функционального стиля тогдашнего constexpr-кода; но это бы фактически создало новый под-язык constexpr C++ функционального стиля.

Самым близким к "остальному C++" вариантом было не заменять качество количеством, а просто стараться поддерживать в constexpr-вычислениях широкое подмножество C++ (в идеале - его весь). Этот вариант был выбран. Это сильно повлияло на дальнейшую историю constexpr.

Поэтому возникла необходимость в мутабельности объектов в рамках constexpr-вычислений. Согласно предложению, объект, созданный в рамках вычисления constexpr-выражения, теперь можно изменять в течение процесса вычисления, до тех пор, пока не закончится процесс вычисления или лайфтайм объекта.

Эти вычисления до сих пор происходят внутри своей "песочницы", ничто снаружи на них не влияет, поэтому, по идее, вычисление constexpr-выражения с одними и теми же аргументами будет давать один и тот же результат (не считая погрешностей в float- и double-вычислениях).

Для лучшего понимания я скопировал кусок кода из предложения:

constexpr int f(int a) {
  int n = a;
  ++n;                  // '++n' is not a constant expression
  return n * a;
}

int k = f(4);           // OK, this is a constant expression.
                        // 'n' in 'f' can be modified because its lifetime
                        // began during the evaluation of the expression.

constexpr int k2 = ++k; // error, not a constant expression, cannot modify
                        // 'k' because its lifetime did not begin within
                        // this expression.

struct X {
  constexpr X() : n(5) {
    n *= 2;             // not a constant expression
  }
  int n;
};
constexpr int g() {
  X x;                  // initialization of 'x' is a constant expression
  return x.n;
}
constexpr int k3 = g(); // OK, this is a constant expression.
                        // 'x.n' can be modified because the lifetime of
                        // 'x' began during the evaluation of 'g()'.

От себя замечу, что теперь компилируется такой код:

constexpr void add(X& x) {
    x.n++;
}
constexpr int g() {
  X x;
  add(x);
  return x.n;
}

Теперь в constexpr-методах работает значимая часть С++, и в методах разрешены сайд-эффекты, локальные в рамках constexpr-вычисления. Вычислитель констант стал сложнее, но все еще справлялся с задачей.

2013: Легендарные const-методы и популярные constexpr-методы

constexpr-методы класса на данный момент автоматически помечаются как const-методы.

В предложении [N3598] constexpr member functions and implicit const обратили внимание, что constexpr-методы класса не обязательно неявно делать const-методами.

Это стало актуальнее с мутабельностью в constexpr-вычислениях; но и до этого мешало использовать один и тот же метод в constexpr и не-constexpr коде:

struct B {
  constexpr B() : a() {}
  constexpr const A &getA() const /*implicit*/ { return a; }
  A &getA() { return a; } // дублирование кода
  A a;
};

Интересно, что предложение давало на выбор три опции, из них был выбран второй:

Статус-кво; минус: дублирование кода
constexpr не будет неявно значить const; минус: ломает ABI (const является частью mangled-имени метода)
Добавить новый квалификатор и писать constexpr A &getA() mutable { return a; }; минус: новый баззворд в конце объявления

2015-2016: Синтаксический сахар для шаблонов

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

template <class T, class... Args> 
enable_if_t<is_constructible_v<T, Args...>, unique_ptr<T>> 
make_unique(Args&&... args) 
{
    return unique_ptr<T>(new T(forward<Args>(args)...));
}  

template <class T, class... Args>  
enable_if_t<!is_constructible_v<T, Args...>, unique_ptr<T>>
make_unique(Args&&... args) 
{
    return unique_ptr<T>(new T{forward<Args>(args)...});
}

Предложение [N4461] Static if resurrected вводит выражение static_if (позаимствовав из языка D), чтобы код стал менее страшным:

template <class T, class... Args> 
unique_ptr<T>
make_unique(Args&&... args) 
{
    static_if (is_constructible_v<T, Args...>) {
        return unique_ptr<T>(new T(forward<Args>(args)...));
    } else {
        return unique_ptr<T>(new T{forward<Args>(args)...});
    }
}

Этот кусок C++ имеет довольно посредственное отношение к constexpr-вычислениям и работает в другой сфере. Но static_if в следующих ревизиях решили переименовать:

constexpr_if (is_constructible_v<T, Args...>) {
  return unique_ptr<T>(new T(forward<Args>(args)...));
} constexpr_else {
  return unique_ptr<T>(new T{forward<Args>(args)...});
}

Затем еще немного...

constexpr if (is_constructible_v<T, Args...>) {
        return unique_ptr<T>(new T(forward<Args>(args)...));
} constexpr_else {
        return unique_ptr<T>(new T{forward<Args>(args)...});
}

И финальный вариант:

if constexpr (is_constructible_v<T, Args...>) {
    return unique_ptr<T>(new T(forward<Args>(args)...));
} else {
    return unique_ptr<T>(new T{forward<Args>(args)...});
}

2015: Constexpr-лямбды

В очень хорошем предложении [N4487] Constexpr Lambda подробно проработали вопрос использования closure type в constexpr-вычислениях (и написали поддержку в форкнутом clang).

Чтобы понять, как возможно иметь constexpr-лямбды, нужно понимать, как они устроены "внутри". Есть статья про историю лямбд, где описано, как прото-лямбды существовали уже в C++03, и в сегодняшних лямбда-выражениях создается похожий класс, скрытый за чертогами компилятора.

Прото-лямбда для [](int x) { std::cout << x << std::endl; }

#include <iostream>
#include <algorithm>
#include <vector>

struct PrintFunctor {
    void operator()(int x) const {
        std::cout << x << std::endl;
    }
};
int main() {
    std::vector<int> v;
    v.push_back(1);
    v.push_back(2);
    std::for_each(v.begin(), v.end(), PrintFunctor());
}

Если все "захваченные" переменные являются литеральными типами, то closure type тоже предложено считать литеральным типом, а operator() пометить constexpr. Работающий пример constexpr-лямбд:

constexpr auto add = [] (int n, int m) {
    auto L = [=] { return n; };
    auto R = [=] { return m; };
    return [=] { return L() + R(); };
};
static_assert(add(3, 4)() == 7, "");

2017-2019: Двойные стандарты

В предложении [P0595] The constexpr Operator рассмотрели возможность "знать" внутри метода, где сейчас выполняется метод - в вычислителе констант или в рантайме. Автор предложил использовать для этого вызов constexpr(), которые будет соответственно возвращать true или false.

constexpr double hard_math_function(double b, int x) {
  if (constexpr() && x >= 0) {
    // медленная формула, более точная (compile-time)
  } else {
    // быстрая формула, менее точная (run-time)
  }
}

Затем оператор был заменен на "магическую" функцию std::is_constant_evaluated() ([P0595R2]) и в таком виде принят в Стандарт С++20.

Если предложение разрабатывается долго, то авторы иногда делают "rebase" (аналогично как в проектах в git/svn), приводя его в соответствие с обновившимся состоянием.

Так и здесь - авторы [P1938] if consteval (про consteval будет позже) обнаружили, что лучше создать новую запись:

if consteval { }
if (std::is_constant_evaluated()) { }
// ^^^ аналогичные записи

Это решение было принято в C++23 - ссылка на голосование.

2017-2019: We need to go deeper

В constexpr-методах во время constexpr-вычислений пока нельзя использовать дебаггер и выводить логи. Предложение [P0596] std::constexpr_trace and std::constexpr_assert рассматривает введение специальных методов для этих целей.

Это предложение было благосклонно принято - ссылка на голосование, но пока еще не доработано.

2017: Злой двойник стандартной библиотеки

На данный момент std::vector, который желательно затащить в compile-time, не может работать в constexpr-вычислениях, в основном из-за недоступности там операторов new/delete.

Идея о допуске операторов new и delete в вычислитель констант на тот момент выглядела слишком амбициозно, поэтому в довольно странном предложении [P0597] std::constexpr_vector рассматривается введение магического std::constexpr_vector<T>.

Он является противоположносью std::vector<T> - может быть создан и изменен только во время constexpr-вычислений.

constexpr constexpr_vector<int> x;  // Okay.
constexpr constexpr_vector<int> y{ 1, 2, 3 };  // Okay.
const constexpr_vector<int> xe;  // Invalid: not constexpr

Как вычислитель констант должен работать с памятью - не описано. @antoshkka и @ZaMaZaN4iK (авторы многих предложений) в [P0639R0] Changing attack vector of the constexpr_vector выявили большое количество минусов подхода, и предложили поменять направление работы в сторону абстрактного магического constexpr allocator, который не дублирует всю стандартную библиотеку.

2017-2019: Constexpr обретает память

В презентации Constexpr ALL the thing! демонстрировался пример constexpr-библиотеки для работы с JSON-объектами (то же самое, но в бумажном виде, есть в [P0810] constexpr in practice):

constexpr auto jsv
    = R"({
          "feature-x-enabled": true,
          "value-of-y": 1729,
          "z-options": {"a": null,
                        "b": "220 and 284",
                        "c": [6, 28, 496]}
         })"_json;
if constexpr (jsv["feature-x-enabled"]) {
    // code for feature x
} else {
    // code when feature x turned off
}

Авторы сильно пострадали от невозможности использовать контейнеры STL, и написали свои велосипедные аналоги std::vector и std::map, имеющие внутри себя std::array (умеющий работать в constexpr).

Предложение [P0784] Standard containers and constexpr исследует возможность ввода STL-контейнеров в constexpr-вычисления.

[Note: Важно знать, что такое аллокатор. STL-контейнеры работают через него с памятью. Какой именно аллокатор - задается через аргумент шаблона. Чтобы войти в тему, можно почитать эту статью. end note]

Что же мешает разрешить STL-контейнеры в constexpr-вычислениях? Есть три проблемы:

Деструкторы не могут быть объявлены constexpr (у constexpr-объектов он обязан быть тривиальным).
Недоступна динамическая аллокация/деаллокация памяти.
Недоступен placement-new для вызова конструктора в аллоцированной памяти.

Первая проблема. С этим разобрались быстро - авторы предложения обсуждали проблему с разработчиками фронтенда MSVC++, GCC, Clang, EDG, и они подтвердили, что ограничение можно легко ослабить. Теперь от литеральных типов можно требовать наличие constexpr-деструктора, а не строго тривиального деструктора.

Вторая проблема. Работа с памятью не очень проста. Как уже упоминалось, вычислитель констант обязан отлавливать undefined behaviour в любом виде и останавливать компиляцию в случае его наличия.

Это значит, что вместе с многими объектами нужно трекать их "метаданные", которые держат руку на пульсе и не дают сломать выполнение программы. Пара примеров таких метаданных:

Информация, какое поле в union-е активно ([P1330]) (примеры undefined behavior: запись в член неактивного поля)
Жесткая связь между указателем или ссылкой и соответствующим ему реальным ранее созданным объектом (примеры undefined behavior: бесконечное множество)

Из-за этого не имеет смысла использовать подобные методы:

void* operator new(std::size_t);

Так как нет никакого обоснования приводить void* к T*. Короче говоря, новая ссылка/указатель либо могут начать указывать на существующий объект, либо быть созданным "одновременно" с ним, и никак иначе.

Поэтому есть две опции работы с памятью, которые будут законны в constexpr-вычислениях:

Простые new- и delete-выражения: int* i = new int(42);
Использование стандартного аллокатора: std::allocator (его подпилили напильником)

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

Контейнеры полагаются на std::allocator_traits, в частности для конструирования - на его метод construct. До предложения он имел вид

template< class T, class... Args >
static void construct( Alloc& a, T* p, Args&&... args ) {
        ::new (static_cast<void*>(p)) T(std::forward<Args>(args)...);
        // ^^^ placement-new, запрещенный в constexpr-вычислениях
}

То есть его нельзя использовать из-за приведения к void* и placement-new (которые в общем виде в constexpr запрещены). А в предложении преобразовался в

template< class T, class... Args >
static constexpr void construct( Alloc& a, T* p, Args&&... args ) {
        std::construct_at(p, std::forward<Args>(args)...);
}

std::construct_at - это метод, который в рантайме работает аналогично старому коду (с приведением к void*), а в constexpr-вычислениях он:

.∧＿∧
( ･ω･｡)つ━☆・*。
⊂　 ノ 　　　・゜+.
しーＪ　　　°。+ *´¨)
　　　　　　　　　.· ´¸.·*´¨) ¸.·*¨)
　　　　　　　　　　(¸.·´ (¸.·'* ☆ Вжуххххх, и просто работает! ☆

То есть компиляторный вычислитель констант обработает его особым образом: по видимости, вызвав конструктор у объекта, который связан с T* p, без бюрократии.

Этого достаточно, чтобы контейнеры стало возможно использовать в constexpr-вычислениях.

На первых порах на аллоцированную память наложили ограничения: она должна быть деаллоцирована в рамках этого же constexpr-вычисления, не выходя за рамки "песочницы".

Такой новый тип аллокации памяти назван transient constexpr allocations. Слово transient приблизительно можно перевести как "мимолетный" или "недолговечный".

В этом предложении еще был кусок про non-transient allocation - дать возможность освобождать не всю аллоцированную память, тогда неосвобожденная память "вываливается" из "песочницы" и конвертируется бы в static storage (другими словами, в секцию .rodata). Но Комитет посчитал эту возможность "слишком хрупкой" ("too brittle") по многим причинам, и пока его не принял.

В остальном, это предложение было принято.

2018: Поймай меня, если сможешь

Предложение [P1002] Try-catch blocks in constexpr functions вводит try-catch блоки в constexpr-вычисления.

Это предложение немного сбивает с толку, так как throw на тот момент был запрещен в constexpr-вычислениях (значит, catch-кусок кода никогда не запускается).

Судя по документу, это ввели, чтобы пометить все методы std::vector как constexpr - в libc++ (реализация STL) в методе vector::insert используется try-catch блок.

2018: Я сказал constexpr!

Из личного опыта знакомо, что двойственность constexpr-методов (могут выполняться и в compile-time, и в run-time) приводит к тому, что вычисления проваливаются в рантайм там, где не ожидаешь: пример кода; и для того, чтобы гарантировать нужный этап, нужно фокусничать: пример кода.

Предложение [P1073] constexpr! functions вводит новое ключевое слово constexpr! для методов, которые должны работать только в compile-time. Эти методы называются immediate-методами.

constexpr! int sqr(int n) {
  return n*n;
}
constexpr int r = sqr(100);  // Okay.
int x = 100;
int r2 = sqr(x);  // Error: Call does not produce a constant.

Если существует возможность того, что в constexpr!-метод могут попасть переменные, неизвестные на этапе компиляции (что норма для constexpr-методов), то программа не скомпилируется:

constexpr! int sqrsqr(int n) {
  return sqr(sqr(n)); // Not a constant-expression at this  point,
}                     // but that's okay.

constexpr int dblsqr(int n) {
  return 2*sqr(n); // Error: Enclosing function is not
}                  // constexpr!.

К constexpr!-методу нельзя взять указатель или ссылку. Бэкенду компилятора вовсе не обязательно (и не нужно) знать про существование таких функций, помещать их в таблицы символов и т.д.

В следующих ревизиях этого предложения ключевое слово constexpr! поменяли на consteval.

Разница между constexpr и consteval налицо - во втором случае нет провала в рантайм: пример с constexpr, пример с consteval.

2018: Слишком радикальный constexpr

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

Предложение [P1235] Implicit constexpr предлагает по умолчанию помечать все методы, имеющие определение, как constexpr. Но можно и запретить выполнять метод в compile-time:

<нет спецификатора> - метод помечается как constexpr, если возможно.
constexpr - работает как сейчас
constexpr(false) - не может быть вызван в compile-time
constexpr(true) - может быть вызван только в compile-time, т.е. аналогично constexpr!/consteval

Это предложение не было принято - ссылка на голосование.

2020: Долговечная constexpr-память

Как уже обсуждалось, после принятия предложения [P0784] Standard containers and constexpr в constexpr-вычислениях стало возможно аллоцировать память, но ее всю необходимо освобождать до окончания constexpr-вычисления: это так называемые transient constexpr allocations.

Таким образом, нельзя создавать top-level constexpr-объекты почти всех STL-контейнеров и многих других классов.

[Note: Под "top-level объектом" я имею в виду результат всего constexpr-вычисления, например:

constexpr TFoo CalcFoo();
constexpr TFoo FooObj = CalcFoo();

Здесь вызов CalcFoo() начинает constexpr-вычисление, а FooObj - его результат и top-level constexpr-объект. end note]

Предложение [P1974] Non-transient constexpr allocation using propconst находит путь к решению проблемы. На мой взгляд, это самое интересное предложение из всех, что я привел в статье (оно заслуживало бы отдельной статьи). Ему дали "зеленый свет" и оно развивается - ссылка на тикет. Здесь я его перескажу в понятной форме.

Что же мешает нам иметь non-transient allocations? Вообще, проблема не в том, чтобы запихать куски памяти в static storage (.bss/.rodata/их аналоги), а в том, чтобы проверить, что вся схема обладает чёткой консистентностью.

Допустим, что мы имеем некий constexpr-объект, конструирование (точнее, "вычисление") которого спровоцировало non-transient allocations. Значит, теоретическое деконструирование этого объекта (т.е. вызов его деструктора) должно освободить всю non-transient память. Если вызов деструктора не освободил бы память, то это плохо (консистентности нет) и надо выдавать ошибку компиляции.

Другими словами, вот что должен делать вычислитель контант:

Увидев запрос на constexpr-вычисление, произвести его.
В результате вычисления получить объект, скрывающий под собой пачку constexpr-переменных литерального типа; и некоторый объём неосвобожденной памяти (non-transient allocations).
Имитировать вызов деструктора у данного объекта (не вызывая его на самом деле), и проверить что этот вызов освободил бы всю non-transient память.
Если все проверки прошли успешно, то консистентность доказана. Non-transient allocations можно двигать в static storage.

Это звучит логично, и допустим, что это всё реализовано. Но тогда получим проблему с подобным кодом с non-transient памятью, которую стандарт не будет запрещать менять, и тогда проверка на вызов деструктора будет бессмысленной:

constexpr unique_ptr<unique_ptr<int>> uui
    = make_unique<unique_ptr<int>>(make_unique<int>());
int main() {
    unique_ptr<int>& ui = *uui;
    ui.reset();
}

[Note: В реальности такой код получил бы отпор от ОС за попытку записи в read-only сегмент RAM, но это физическая константность. А в коде должна быть логическая константность. end note]

Пометка constexpr у объектов влечет за собой пометку их как const, и все их члены также становятся const.

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

У указательных типов есть два ортогональных параметра константности:

Можно ли начать указывать на другой объект?
Можно ли изменять объект, на который указывается?

И получается 4 варианта с разными свойствами (OK - строка компилируется, FAIL - нет):

int dummy = 13;

int * test1{nullptr};
test1 = &dummy; // OK
*test1 = dummy; // OK

int const * test2{nullptr};
test2 = &dummy; // OK
*test2 = dummy; // FAIL

int * const test3{nullptr};
test3 = &dummy; // FAIL
*test3 = dummy; // OK

int const * const test4{nullptr};
test4 = &dummy; // FAIL
*test4 = dummy; // FAIL

"Обычный" const приводит к третьему варианту, но для constexpr необходим четвертый! То есть необходим так называемый deep-const.

Предложение на базе пары других, более старых предложений, предлагает для этого ввести новый cv-qualifier propconst (propagating const).

Этот квалификатор будет использоваться с указательными/ссылочными типами:

T propconst*
T propconst&

И в зависимости от типа T компилятор будет либо конвертировать это слово в const, либо удалять его. Первый случай - если T константный, второй - если нет:

int propconst * ---> int *
int propconst * const ---> int const * const

В предложении приведена таблица конвертации propconst в разных случаях.

Таким образом, constexpr-объекты могли бы обрести полную логическую константность (deep-const):

constexpr unique_ptr<unique_ptr<int propconst> propconst> uui
    = make_unique<unique_ptr<int propconst> propconst>(make_unique<int propconst>());

int main() {
    unique_ptr<int propconst>& ui = *uui;
    ui.reset();
    // ^^^ FAIL

    const unique_ptr<int propconst>& ui = *uui;
    // ui.reset();
    // ^^^ OK
}

// P.S. Такая запись еще не принята Комитетом, я надеюсь, сделают лучше

2021: Constexpr-классы

С появлением полностью constexpr-классов, включая std::vector, std::string, std::unique_ptr, в которых все методы методы помечены как constexpr, возникает желание сказать "пометь все методы класса как constexpr".

Это делает предложение [P2350] constexpr class:

class SomeType {
public:
  constexpr bool empty() const { /* */ }
  constexpr auto size() const { /* */ }
  constexpr void clear() { /* */ }
  // ...
};
// ^^^ ДО

class SomeType constexpr {
public:
  bool empty() const { /* */ }
  auto size() const { /* */ }
  void clear() { /* */ }
  // ...
};
// ^^^ ПОСЛЕ

С этим предложением связана интересная история - еще не зная о его существовании, я вкинул в stdcpp.ru идею предложить такую же штуку: ссылка на тикет (что сейчас не нужно).

Многие почти одинаковые предложения в Стандарт могут появиться почти одновременно. Это говорит в пользу теории множественных открытий: идеи "витают в воздухе", и в принципе не так важно, кто их открывает - если комьюнити достаточно большое, то естественная эволюция происходит.

2019-∞: Интерпретатор констант в компиляторе

constexpr-вычисления могут быть очень медленными, потому что вычислитель констант на синтаксическом дереве развивался итеративным образом (начиная с constant folding) и сейчас делает множество лишних вещей, которые можно было делать более эффективно.

С 2019 года в Clang разрабывается ConstantInterpeter, который в перспективе может полностью заменить вычислитель констант на синтаксическом дереве. Он довольно интересен и заслуживал бы отдельной статьи.

Его идея состоит в том, что на основе синтаксического дерева можно сгенерировать "байткод", который затем выполнить на интерпретаторе. Интерпретатор поддерживает в себе стек, фреймы вызовов, модель памяти (с метаданными, о которых говорилось ранее).

Документация для ConstantInterpeter хорошая, и также много интересного есть в видео-выступлении создателя интерпретатора на конференции LLVM-разработчиков.

Что можно еще посмотреть?

Для того, чтобы больше расширить понимание, можно посмотреть замечательные выступления от экспертов. В каждом выступлении авторы выходят за рамки рассказа про constexpr: это может быть построение constexpr-библиотеки; рассказ про использование constexpr в будущем reflexpr; или про устройство вычислителя констант и интерпретатора констант.

constexpr ALL the things!, Ben Deane & Jason Turner, C++Now 2017. Уже немного устарело, но может быть интересно, про построение constexpr-библиотеки.
Compile-time programming and reflection in C++20 and beyond, Louis Dionne, CppCon 2018. Много внимания уделяется будущей рефлексии в C++.
Полезный constexpr, Антон Полухин (@antoshkka), C++ CoreHard Autumn 2018. Есть про компиляторы, рефлексию и метаклассы.
The clang constexpr interpreter, Nandor Licker, 2019 LLVM Developers’ Meeting. Рокет саенс и интерпретатор кода для constexpr.

Я хотел бы, чтобы здесь было выступление про киллер-фичу (по моему мнению) [P1040] std::embed, которая отлично работала бы в тандеме с constexpr. Но, судя по тикету, его планируют реализовать в C++-каком-то.

Adblock test (Why?)

Игра в собственные

Имеем набор данных в виде совокупности квадратных матриц, которые используются - вместе с известным выходом - в качестве тренировочного набора для нейронной сети. Можно ли обучить нейронную сеть, используя только собственные значения матриц? Во избежание проблем с комплексными значениями, упор делаем на симметричные матрицы. Для иллюстрации используем набор данных MNIST. Понятно, что невозможно восстановить матрицу по ее собственными значениям - для этого понадобится еще кое-что, о чем мы поговорим далее. Поэтому трудно ожидать некоего прорыва на данном пути, хотя известно, что можно говорить о чем угодно, строить грандиозные планы, пока не пришло время платить. О деньгах мы здесь не говорим, просто задаем глупый вопрос, на который постараемся получить осмысленный ответ, тем более что в процессе познания расширим свои научные горизонты. Например, сначала мы познакомимся с тем, как находить собственные векторы и собственные значения (eigenvalues and eigenvectors) для заданной квадратной матрицы, затем плавно выкатим на эрмитовы и унитарные матрицы. Все иллюстративные примеры сопровождаются простыми кодами. Далее возьмем MNIST , преобразуем в набор собственных значений симметричных матриц и используем молоток от Keras. Как говорят в Японии: “Торчащий гвоздь забивают”. Закроем глаза и начнем бить, а на результат можно и не смотреть: получится как всегда. Сразу скажу, что изложение будет проведено как можно ближе к тому, как я это дело понимаю для себя, не обращаясь к строгому обоснованию, которое обычно не используется в повседневной жизни. Иными словами, что понятно одному глупцу, понятно и другому. Все мы невежественны, но, к счастью, не в одинаковой степени. С другой стороны, предполагаю, что многие, хоть и в гимназиях не обучались, но имеют представление - по своему опыту обучения, - что значит впихнуть невпихуемое.

Собственно, собственные

Итак, пусть у нас задана произвольная квадратная матрица с компонентами $A_{i,j}$ . Мы ставим запятую между индексами, имея в виду оформление массивов в большинстве языков программирования; обычно так не делается, но мы делаем так, как нам хочется, да и выглядит все более разреженным и удобочитаемым. Мы можем умножить матрицу на произвольный вектор v_i , руководствуясь правилом умножения матриц (строка на столбец), $A_{i,j}v_j$ . По повторяющимся индексам производится суммирование; все получится, если размерность первой (нумерация с нуля) оси матрицы совпадает с размерностью вектора. Если будет желание, когда я что-то упускаю из виду, можно заглянуть в мои заметки . Вдруг приключилось так, что $A_{i,j}v_j=\lambda v_i\quad\Rightarrow\quad A\cdot v=\lambda v$ , тогда говорят: $\lambda$ - собственное значение, - собственный вектор. Опять же, возвращаясь к глубокомысленной дискуссии о запятых в индексах, здесь мы используем точку, когда переходим от компонент к матрицам, имея в виду замечательную функцию dot() в превосходном NumPy. Теперь проведем ряд очевидных манипуляций

$A_{i,j}v_j=\lambda v_i\Rightarrow A_{i,j}v_j=\lambda \delta_{i,j} v_j\Rightarrow (A_{i,j}-\lambda \delta_{i,j}) v_j=0\Rightarrow (A-\lambda I)\cdot v=0\tag{1},$

где $\delta_{i,j}$ - символ Кронекера, - единичная матрица, 0 - нулевой вектор. Условие нетривиальности решения, $v\neq 0$ , для собственного вектора приводит к характеристическому уравнению $det(A-\lambda I)=0$ . Если размер матрицы равен n, получим алгебраическое уравнение n-й степени, которое имеет n корней (действительных или комплексных) - собственных значений (eigenvalues) $\lambda_i$ . Каждому собственному значению соответствует собственный вектор (eigenvectors) $v_0, v_i, ...v_{n-1}$ (например, $A\cdot v_3=\lambda_3 v_3$ ). Совокупность собственных векторов образует матрицу с компонентами $v_{i,j}$ , где нулевая ось - компоненты вектора, а первая - нумерация собственных векторов. Теперь уравнение на собственные значения выглядит следующим образом: $A_{i,j}v_{j,k}=\lambda_k v_{i, k}$ . Раз есть вектор, есть его длина, которая равна сумме квадратов компонент (квадрат длины). Например, для квадрата длины третьего вектора (нумерация с нуля) имеем

$v_{i,2}v_{i,2}=(V.T)_{2,i}v_{i,2}\quad\Rightarrow\quad (V.T \cdot V)_{2,2},\tag{2}$

где, вспоминая правила умножения матриц ( $C_{i,k}=A_{i, j}B_{j, k}$ , см. заметки ), мы воспользовались тождеством $v_{i,2}=(V.T)_{2,i}$ , в котором V.T - транспонированная матрица (используем обозначения NumPy). Здесь нет проблем, но это до тех пор, пока мы работаем с действительными векторами. В комплексной области квадрат длины вектора равен скалярному произведению вектора на комплексно сопряженный. В этом случае потребуется следующая модификация соотношения (2):

$v_{i,2}^*v_{i,2}=((V.T)^*_{2,i}v_{i,2}\quad\Rightarrow\quad (V.H \cdot V)_{2,2},\tag{3}$

где звездочка обозначает обычное комплексное сопряжение. Следует заметить, что (V.T)^*=V.H - эрмитово-сопряжённая матрица. Длина вектора - число, отсюда открывается возможность выделить собственное значение путем умножения на собственный вектор. В результате для выделенного собственного значения получим

$A_{i,j}v_{j,2}=\lambda_2 v_{i, 2}\quad\Rightarrow\quad v_{i,2}^*A_{i,j}v_{j,2}=\lambda_2 v_{i,2}^*v_{i,2}\quad\Rightarrow\quad \lambda_2 =\frac{v_{i,2}^*A_{i,j}v_{j,2}}{v_{i,2}^*v_{i,2}}.\tag{4}$

Все это выглядит не очень хорошо, поскольку произведение разных собственных векторов, например $v_{i,2}^*v_{i,3}$ , может быть не равно нулю. Это настораживает, поскольку мы привыкли иметь дело с единичными векторами, каждый из которых ортогонален другому. Вспомните хотя бы единичные векторы по трем осям в обычном пространстве. Но проблемы отложим на потом, а пока создадим видимость работы - побалуемся с кодами.

import numpy as np
from numpy import matrix
from numpy.linalg import eig 
# из субмодуля (submodule linalg) линейной алгебры берем метод
#для вычисления собственных значений и векторов 

# Возьмем матрицу, которую можно проверить вручную
A = matrix([[3, -0.65], [5, 1]])

# calculate eigendecomposition
values, V = eig(A)
print(values) # [2.+1.5j 2.-1.5j]
print(V) # [[0.18814417+0.28221626j 0.18814417-0.28221626j]
         #  [0.94072087+0.j         0.94072087-0.j        ]]

Действительно, характеристическое уравнение $(3-\lambda)(1-\lambda)+3.25=0$ имеет решения $\lambda_{+,-}=2\pm 1.5j$ , где - мнимая единица. Столбцы матрицы - собственные векторы v_0 и v_1 .

# Первый собственный вектор - все элементы нулевого столбца и т.д.
v0=V[:,0]
v1=V[:,1]
v0.shape, v1.shape # ((2, 1), (2, 1))
# К счастью, они единичные
np.dot(v0.H,v0) # matrix([[1.+0.j]])
np.dot(v1.H,v1)[0,0] # (0.9999999999999999+0j). Ну, почти единица
# Но эти векторы не ортогональны
np.dot(v0.H,v1) # matrix([[0.84070796-0.10619469j]])

Эрмитовы матрицы

Мы научились находить eigendecomposition для произвольной квадратной матрицы. Пока не будем говорить о вырожденных матрицах , которых, по возможности, будем в дальнейшем избегать. Итак, в чем состоит загвоздка. С собственными значениями - все в порядке, а вот собственные векторы изрядно подгуляли. Дело в том, что они не ортогональны друг другу. Посмотрим, какие условия на матрицу нам потребуется наложить, чтобы добиться ортогональности, заодно откроем для себя ряд новых возможностей в деле познания собственных способностей.

Пусть из всего набора мы произвольно выделили два участника $A_{i,j}v_{j,k}=\lambda_k v_{i, k}$ и $A_{i,j}v_{j,l}=\lambda_l v_{i, l}$ , где, напомним, второй индекс нумерует собственные вектора. Берем первое соотношение и умножаем его на $v_{i,l}^*$ , а второе - на $v_{i,k}^*$ , в результате имеем

$v_{i,l}^*A_{i,j}v_{j,k}=\lambda_k v_{i,l}^*v_{i, k}\quad\quad v_{i,k}^*A_{i,j}v_{j,l}=\lambda_l v_{i,k}^*v_{i, l}.\tag{5}$

Далее выполним комплексное сопряжение второго уравнения в (5)

$v_{i,k}A_{i,j}^*v_{j,l}^*=\lambda_l^* v_{i,k}v_{i, l}^*\quad\Rightarrow\quad v_{j,l}^*A_{i,j}^*v_{i,k}=\lambda_l^* v_{i,l}^*v_{i, k}.$

Поскольку индексы в последнем уравнении “немые” (по ним производится суммирование) , их можно назвать как угодно. Таким образом, после замены $i\Leftrightarrow j$ получим

$v_{i,l}^*A_{j,i}^*v_{j,k}=\lambda_l^* v_{i,l}^*v_{i, k}\quad\Rightarrow\quad v_{i,l}^*(A.H)_{i,j}v_{j,k}=\lambda_l^* v_{i,l}^*v_{i, k},\tag{6}$

где A.H - эрмитово-сопряжённая матрица, элементы которой определяются как $(A.H)_{i,j}=A_{j,i}^*$ . Теперь из первого уравнения (5) вычитаем (6). В результате получим следующее замечательное выражение:

$v_{i,l}^*[A-A.H]_{i,j}v_{j,k}=(\lambda_k-\lambda_l^*) v_{i,l}^*v_{i, k}.\tag{7}$

Пока мы не накладывали никаких ограничений на матрицу. Первое, что приходит на ум: положим A=A.H , определяя тем самым эрмитову матрицу Следует отметить: все что мы наблюдаем в окружающем нас мире связано, так или иначе, с собственными значениями эрмитовых матриц. Шарль Эрмит - вот кто Создатель Вселенных (см. превосходную биографию, Ожигова Е.П. Шарль Эрмит. - Л.:Наука, 1982.). Да, было время, когда математики занимались реальным миром, не то, что нынешнее племя.

Итак, для эрмитовых матриц получаем выражение $(\lambda_k-\lambda_l^*) v_{i,l}^*v_{i, k}=0$ , которое связано с двумя важнейшими событиями для каждого, кто изучал квантовую механику. Во-первых, собственные значения эрмитовых матриц являются действительными величинами , поскольку при l=k из условия $v_{i,k}^*v_{i, k}>0$ имеем $\lambda_k=\lambda_k^*$ . Во- вторых, если $l\neq k$ и $\lambda_k\neq\lambda_l$ , получаем свойство ортогональности собственных векторов $v_{i,l}^*v_{i, k}=0$ . Для эрмитовых матриц собственные значения - действительные величины, а собственные векторы - ортогональны друг другу, вернее, ортонормированы, поскольку удобно включать условие нормировки. Таким образом, мы стремились к аналогии с обычным пространством, с единичными векторами по ортогональным осям, мы это и получили.

Теперь соорудим эрмитову матрицу. Для этого сгенерируем пару случайных матриц и , затем сложим их, умножив одну из матриц на мнимую единицу , A=X+jY . Затем возьмем полученную матрицу и сложим ее с эрмитово-сопряженной. Действительно, $B=A+A.H\Rightarrow B.H=(A+A.H).H=A.H+A=B$ , так что - эрмитова матрица.

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

X=np.random.randint(0, 256, (28,28))# Вспоминаем молодость с MNIST
Y=np.random.randint(0, 256, (28,28))
A=matrix(X+1j*Y)
B=A+A.H
# выводим собственные значения и собственные вектора (в виде матрицы)
values, V = eig(B)
values=np.real(values)# только действительная часть
values

Собственные значения (values) получились с едва заметной мнимой частью, которую мы отрезаем по идеологическим соображениям. Для этого мы оставили только действительную часть. Напомним, что

$v_{i,l}^*v_{i, k}=(V.H)_{l,i}v_{i, k}=\delta_{l,k}\quad\Rightarrow\quad V.H\cdot V=E,\tag{8}$

где - единичная матрица (по диагонали единицы ). Проверим в коде

np.dot(V.H, V) # почти единичная матрица с формой (28,28)

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

Симметричный MNIST

Итак, в предыдущем разделе мы говорили о произвольных о эрмитовых матрицах. Если же мы имеем дело с действительными матрицами, комплексное сопряжение выпадает, поэтому эрмитовы - это просто симметричные матрицы, A=A.T .

В этом разделе мы возьмем набор данных MNIST и попробуем поработать с каждым элементом как с матрицей. Каждое рукописное изображение цифр (от 0 до 9) - матрица с формой (28,28), каждый элемент которой - число от 0 до 255 (8-битная шкала серого). Загрузим и посмотрим (рекомендую начать новый блокнот).

import keras
from keras.datasets import mnist
keras.__version__ # '2.4.3'
(x,y),(tx,ty)=mnist.load_data()

База данных MNIST состоит тренировочного (x,y) и тестового (tx,ty) наборов, где - набор рукописных изображений цифр(матриц), которые подаем на вход, а - набор правильных ответов (числа от 0 до 9); тестовый набор выделен для удобства. Теперь посмотрим, сколько их там, и что они из себя представляют. Для этого подключим “рисовалку”.

import matplotlib.pyplot as plt
# Сколько их?
x.shape[0], tx.shape[0] # (60000, 10000)
# Возьмем произвольный элемент
y[12345] # 3
# Посмотрим на каракули
plt.imshow(x[12345], cmap=plt.get_cmap('gray'))
plt.show()

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

x[12345]

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

import numpy as np
from numpy.linalg import eig 
np.linalg.det(x[12345])# 0.0 - детерминант

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

values, V = eig(x[12345])
values # array([   0.          +0.j        ,    0.          +0.j        ,
                1474.46266075  +0.j        ,  838.19932232  +0.j        ,
               -400.18280429+384.94168117j, -400.18280429-384.94168117j,
                297.89145581  +0.j        ,   90.51616055+144.4741198j ,
                 90.51616055-144.4741198j ,  -57.05178416 +68.44040035j,
                -57.05178416 -68.44040035j,    3.05833828 +80.86506799j,
                  3.05833828 -80.86506799j,  -32.62937655 +32.08539788j,
                -32.62937655 -32.08539788j,    4.36245913 +51.20712637j,
                  4.36245913 -51.20712637j,  -16.0035565   +0.j        ,
                  6.13140522  +0.j        ,   59.17272648  +0.j        ,
                  0.          +0.j        ,    0.          +0.j        ,
                  0.          +0.j        ,    0.          +0.j        ,
                  0.          +0.j        ,    0.          +0.j        ,
                  0.          +0.j        ,    0.          +0.j        ])

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

Итак, хочется поскорее избавиться от комплексных чисел. Мы знаем, что эрмитовы матрицы обладают действительными собственными значениями. Поскольку матрица изображения полностью состоит из действительных чисел, то эрмитова матрица - просто симметричная матрица. Иными словами, нам потребуется соорудить симметричную матрицу, но перед этим избавиться от множества нулей. Для этого “закрасим” пикселями, выбирая их значения случайным образом.

from numpy import matrix
a=matrix(x[12345]+np.random.randint(0, 5, (28,28)))
np.linalg.det(a) # У меня -3.7667039553765456e+43, 
                 # у вас, конечно, получится другое значение

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

Итак, перейдем к симметричной матрице и посмотрим, что там нарисовано.

a=a+a.H # для действительных матриц a.H=a.T
plt.imshow(a, cmap=plt.get_cmap('gray'))
plt.show()

Кто скажет, что это не тройка, пусть первый бросит в меня камень! Теперь это число подготовлено для манипуляций в нейронной сети. Матрица симметричная, невырожденная, имеет замечательный набор собственных значений.

values, V = eig(a)
values # array([ 4.06238155e+03,  1.87437218e+03, -1.59889723e+03, -8.51935039e+02,
                5.58412899e+02, -3.53262264e+02, -2.36036195e+02,  2.29960516e+02,
                2.08312088e+02, -1.49974962e+02,  1.40811760e+02, -9.71188401e+01,
                8.27836555e+01, -6.51480016e+01, -4.65323368e+01,  3.90299000e+01,
                3.06552503e+01,  2.06100661e+01,  1.73972579e+01, -1.16491885e+01,
               -8.13295407e+00,  1.24942529e+01, -4.67339958e+00,  9.33795587e+00,
                6.67620627e+00,  2.96767577e+00,  1.13851798e+00,  1.86752842e-02])

Если с собственными значениями все в порядке, то собственные векторы малость подгуляли: мало того, что они комплексные (как обычно и бывает), так восстановление первоначальной матрицы по собственным векторам и собственным значениям может привести к комплексной матрице. Python не знает, что это не правильно, поэтому и вытянет наружу комплексные остатки. Само собой, можно от них избавиться с помощью метода np.real(), оставляя только действительную часть. Но пока попробуем в лоб, проверим прозорливость создателей NumPy на примере восстановления матрицы по собственным значениям и собственным векторам. Чтобы подготовить почву, необходимо вернуться к уравнениям (8), которые справедливы и для нашей матрицы , в которой столбцы - собственные векторы матрицы (см. выше). Таким образом, $а_{i,j}v_{j,k}=\lambda_k v_{i, k}$ . Следует напомнить, что, несмотря на повторяющийся индекс , в правой части нет суммирования (берем фиксированное ). Далее имеем

$а_{i,j}v_{j,k}=v_{i, l}\delta_{l,k}\lambda_k \quad\Rightarrow\quad a\cdot V=V\cdot e(\lambda),$

где $e(\lambda)$ - диагональная матрица собственных значений (по диагонали - собственные значения). Умножая последнее уравнение с правой стороны на V.H , получим

$a\cdot V\cdot V.H=V\cdot e(\lambda)\cdot V.H \quad\Rightarrow\quad a=V\cdot e(\lambda)\cdot V.H,$

в котором использовали (8) . Теперь проверим в коде

e=np.diag(values) # сооружаем диагональную матрицу из собственных значений
b=np.dot(V,np.dot(e,V.H)) # восстанавливаем исходную матрицу по собственным
                          # значениям и векторам; ради приличия ввели новое 
                          # имя матрицы   
plt.imshow(b, cmap=plt.get_cmap('gray'))
plt.show()

Нет слов, все работает! Итак, наши подозрения были беспочвенны, но проверять все равно надо. Теперь мы знаем, что по собственным значениям и векторам можно четко восстановить исходную матрицу. Итак, если используем обратимые операции, царапая формулы на листке бумаге, лучше проверить, как с этим обстоит дело в коде.

Эксперимент

Итак, теперь мы умеем манипулировать симметричными матрицами, используя в качестве полигона базу данных MNIST. Теперь мы попытаемся представить набор данных как симметричные матрицы, затем вычислить собственные значения, на основе которых построить нейронную сеть по примеру обращения с обычным набором. Сразу скажу, что надежд на удачный исход мало. Действительно, как мы выяснили ранее, для того чтобы восстановить матрицу собственных значений недостаточно, для этого необходим набор собственных векторов, объединенных в матрицу , которая сама по себе имеет размер исходной матрицы. Так что, с первого взгляда, мы ничего не выигрываем. Тем не менее, есть надежда, что распределение и свойства собственных значений обладают необходимыми свойствами, чтобы их классифицировать по классам. Иными словами, собственные значения несут информацию о классах рукописных цифр. Распределением собственных значений мы займемся в следующей публикации, а сейчас будем использовать наивный подход, а именно : 1) представим набор MNIST в виде симметричных матриц (28*28); 2) для каждого экземпляра вычислим собственные значения (28 штук); 3) используем Keras. Первые шаги мы уже сделали в предыдущем разделе, а над последним особо заморачиваться не будем, поскольку мы просто ставим эксперименты.

Начинаем новый блокнот:

import numpy as np
from numpy.linalg import eig
import keras
from keras.datasets import mnist

(x,y),(tx,ty)=mnist.load_data()

Теперь напишем функцию, которая будет “закрашивать” матрицы из набора данных. Отмечу, что использую доморощенные функции, которые предпочитаю писать самостоятельно, чем искать на мутных форумах.

def paintSym(z):
    n=z.shape[0]    
    x=np.empty([n,28,28])
    for i in range(n):
        x[i]= z[i]+np.random.randint(0, 5, (28,28))
    return x

Используем эту функцию, затем масштабируем (это не обязательно).

x=paintSym(x)
x=x.astype('float32') / 260

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

def eigSym(z):
    n=z.shape[0]    
    x=np.empty([n,28])
    for i in range(n):
        val, vec = eig((z[i]+z[i].T)/2)
        x[i]= val
    return x

Включаем и смотрим.

x=eigSym(x)
x[12345] # на всякий случай
# array([ 7.82546663e+00,  3.60277987e+00, -3.06879640e+00, -1.62636960e+00,
        1.07726467e+00, -6.77050233e-01, -4.54729229e-01,  4.45481062e-01,
        4.00029063e-01, -2.85397410e-01,  2.71449655e-01, -1.94830164e-01,
        1.41906023e-01, -1.15398742e-01, -8.69008303e-02,  7.79602975e-02,
        5.81430644e-02, -3.25826630e-02,  3.82260419e-02,  3.11152730e-02,
        1.78121049e-02, -1.12444209e-02, -6.83658011e-03,  8.95012729e-03,
        8.10696371e-03,  3.64418700e-03, -2.08872650e-03,  4.40481490e-05])
x.shape # (60000, 28)

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

from keras import models
from keras import layers

network = models.Sequential()
network.add(layers.Dense(512, activation='relu', input_shape=(28,)))
network.add(layers.Dense(10, activation='softmax'))

network.compile(optimizer='rmsprop',
                loss='categorical_crossentropy',
                metrics=['accuracy'])

from keras.utils import to_categorical
y = to_categorical(y)

network.fit(x, y, epochs=100, batch_size=128)

За основу мы берем код из Глубокое обучение на Python. Если осмысленно побалуетесь с цифрами, что-то получите. Точность, которая едва-едва заваливает за 0,5 (за 100 эпох!!!), подкачала по сравнению с оригиналом (accuracy: 0.9897 за Epoch 5), когда используется полная матрица изображения, но хотя бы мы не множили сущности. Для обсуждения не требуется глубокомысленных обобщений, просто, как я наивно думал, если зайца долго бить по голове, он научится спички зажигать. Но оказалось, что этот принцип в глубоком обучении не всегда срабатывает. Возиться дальше - интересно, но бессмысленно. Самое главное, по пути к этому мутному результату можно много чему научиться.

Пару слов напоследок

Итак, пришло время подводить итоги. Понятно, что говорить пока не о чем. Мы попытались на собственных значениях матриц построить работающую схему. Не получилось. Но известно, что отрицательный результат - куда более значимое событие, чем положительный, поскольку при этом у нас открываются новые перспективы. Если бы у нас все получилось, настало бы время для рутинной работы, а так - продолжим идти дальше. Правда и ложь - два взгляда на одну и ту же реальность. Но мы, по крайней мере, не лгали. Просто посмотрели с другого ракурса, заодно открыли для себя новые возможности. К примеру, почему бы не найти преобразование, позволяющее связать все изображения одного класса без использования внешнего наблюдателя для подтверждения события? Иными словами, из одной цифры получить почти полную совокупность из данного набора данных; из “тройки” почти все “тройки”, и т.д.. Сразу приходит на ум конформное отображение двумерных поверхностей: локальные вращения и растяжения с сохранением углов между кривыми, а значит с сохранением формы бесконечно малых фигур. Заодно мы могли бы проследить за трансформацией собственных значений, дабы по их расположению выйти на классификацию по классам. И не обязательно использовать действительные собственные значения. Конечно, переход в комплексную область потребует дополнительных усилий, но там уже существует путь в один путь, многое успешно отработано.

Что касается реальных, на настоящий момент, дел, то это классификация собственных значений по примеру распределения уровней сложных систем. Опять же MNIST животворящий. Но это тема следующей публикации.