У меня возникла идея, как можно расширить синтаксис C++ операцией скалярного произведения. Если кратко, то произведение двух матриц в новых обозначениях будет выглядеть так:
C[>i][>j] = A[i][>k] * B[>k][j];
Насколько мне известно, сочетания операторов [> и [< вроде бы нигде не используются. Их можно применить для декларации индексов, которые существуют только в пределах данного выражения. Сочетание [> используется для декларации индекса, который пробегает от начала до конца массива в прямом направлении, а сочетание [< для декларации индекса, который пробегает в обратном направлении. Для повторяющихся индексов в произведении подразумевается суммирование - они аналогичны немым индексам в тензорных обозначениях.
Разберём на примерах, как это будет работать.
Скалярное произведение двух векторов
float a[3], b[3];
// ...
float result = a[>i] * b[>i];
Здесь мы сначала декларируем конструкцией [>i] индекс i, который будет проходить значения 0,1, 2. После знака произведения повторная конструкция [>i] означает, что этот индекс - немой, и по нему будет идти суммирование.
Данный код эквивалентен следующему:
float a[3], b[3];
// ...
float result = 0;
for(int i = 0; i < 3; i++)
result += a[i] * b[i];
Так как в конструкции a[>i] * b[>i] очевидно, по какому индексу производится суммирование, то допустима и сокращённая запись скалярного произведения, без явного указания имени индекса:
float a[3], b[3];
// ...
float result = a[>] * b[>];
Явное задание размера массива
Так как компилятор не всегда может знать размер динамического массива, то этот размер при декларации индекса придётся указать явно. Это делается при помощи конструкции
[> index: expr]
где результатом выражения expr является размер массива:
float *a = new float[10];
float *b = new float[10];
// ...
float result = a[>i: 10] * b[>i];
Здесь при первой декларации индекса i мы указали, что размер массива равен 10. При последующих декларациях индекса с тем же именем в том же выражении, размер указывать не обязательно.
Данный код эквивалентен следующему:
float *a = new float[10];
float *b = new float[10];
// ...
float result = 0;
for(int i = 0; i < 10; i++)
result += a[i] * b[i];
Поэлементное присвоение вектору выражения
float a[10];
a[>i] = i * i;
Здесь мы декларировали свободный индекс [>i], который говорит компилятору, что мы инициализируем каждый элемент вектора. В данном случае мы присваиваем каждому элементу вектора квадрат его индекса.
Данный код эквивалентен следующему:
float a[10];
for(int i = 0; i < 10; i++)
a[i] = i * i;
Если мы инициализируем вектор выражением, которое от индекса не зависит, то имя индекса можно не указывать:
float a[10];
a[>] = 0;
Поэлементное произведение векторов
float a[3], b[3];
// ...
float result[3];
result[>i] = a[i] * b[i];
В правой части последнего выражения индекс не является немым, и обрабатывается, как обычный индекс массива.
Данный код эквивалентен следующему:
float a[3], b[3];
// ...
float result[3];
for(int i = 0; i < 3; i++)
result[i] = a[i] * b[i];
Реверсия
Формально, выражение присваивается не сразу, а через промежуточный временный объект. Это означает, что пока все итерации не будут завершены - предыдущие значения массивов в правой части выражения не пострадают. Такой подход позволяет простым способом записать переворот вектора (реверсию):
float a[10];
// ...
a[>i] = a[<i];
Здесь в правой части через [<i] мы объявили, что индекс i отсчитывается в обратном направлении: для последнего элемента он будет = 0, для предпоследнего = 1, и так далее.
Данный код эквивалентен следующему:
float a[10];
// ...
for(int i = 0; i < 5; i++)
{
float t = a[i];
a[i] = a[9 - i];
a[9 - i] = t;
}
При реверсии допустимо не указывать явно имя индекса:
float a[10];
// ...
a[>] = a[<];
Умножение матриц
float A[3][3], B[3][3];
// ...
float C[3][3];
C[>i][>j] = A[i][>k] * B[>k][j];
Здесь мы с помощью C[>i][>j] объявляем, что будем инициализировать все элементы матрицы C, затем вычисляем произведение матриц классическим образом: умножая строку матрицы A на столбец матрицы B.
Данный код эквивалентен следующему:
float A[3][3], B[3][3];
// ...
float C[3][3];
for(int i = 0; i < 3; i++)
for(int j = 0; j < 3; j++)
{
float c = 0;
for(int k = 0; k < 3; k++)
c = A[i][k] * B[k][j];
C[i][j] = c;
}
При умножении зачастую выгодно предварительно транспонировать матрицу правого сомножителя. Это тоже просто записывается:
float A[3][3], B[3][3];
// ...
float Bt[3][3];
Bt[>i][>j] = B[j][i]; // транспонирование
float C[3][3];
C[>i][>j] = A[i][>k] * Bt[j][>k]; // эквивалентно A[i][>k] * B[>k][j];
Векторное произведение
Мы можем умножать более двух объектов.Это позволяет записать векторное произведение:
float e[3][3][3]; // псевдотензор Леви-Чивиты
float a[3], b[3];
// ...
e[>][>][>] = 0;
e[0][1][2] = 1;
e[1][0][2] = -1;
e[1][2][0] = 1;
e[2][1][0] = -1;
e[2][0][1] = 1;
e[0][2][1] = -1;
float result[3];
result[>i] = e[i][>j][>k] * a[>j] * b[>k];
Последнее выражение эквивалентно следующему коду:
for(int i = 0; i < 3; i++)
{
float r = 0;
for(int j = 0; j < 3; j++)
for(int k = 0; k < 3; k++)
r += e[i][j][k] * a[j] * b[k];
result[i] = r;
}
Псевдотензор Леви-Чивиты можно сделать объектом стандартной библиотеки. Тогда компилятор сможет обнаруживать его использование, что упростит оптимизацию кода.
Пример из трёхмерной графики
Рассмотрим пусть и не оптимизированный, зато наглядный пример. Предположим, что некий трёхмерный объект задаётся массивом координат вершин, причём вершина, кроме координат, хранит и другие данные. Сначала мы масштабируем трёхмерный объект, причём в каждом направлении осей на свой коэффициент, затем умножаем на матрицу поворота, затем смещаем на какой-то вектор. Результатом получаем новый массив вершин.
// Вершина
struct Vertex
{
float coor[3]; // координаты вершины
// ...
};
int vertexCount; // количество вершин в объекте
float scale[3]; // вектор масштабирующих коэффициентов
float rotation[3][3]; // матрица поворота
float offset[3]; // вектор смещения
// ...
// Создаём исходный объект
Vertex *srcObject = new Vertex[vertexCount];
// ...
// Создаём объект назначения
Vertex *destObject = new Vertex[vertexCount];
destObject[>vertex: vertexCount].coor[>x] =
srcObject[vertex].coor[>i] * scale[i] * rotation[>i][x] + offset[x];
Последнее выражение эквивалентно следующему коду:
for(int vertex = 0; vertex < vertexCount; vertex++)
for(int x = 0; x < 3; x++)
{
float coor = 0;
for(int i = 0; i < 3; i++)
coor += srcObject[vertex].coor[i] * scale[i] * rotation[i][x];
destObject[vertex].coor[x] = coor + offset[x];
}
P.S. Так как в комментариях к этой статье слишком часто писали, что достаточно лишь перегрузить операторы, или ввести функции, которые будут умножать, то я оставлю здесь ссылку на статью о проблемах безындексной записи умножения тензоров: Я не люблю NumPy