Pull to refresh

Возможное расширение языка C++ операцией скалярного произведения

Level of difficultyEasy
Reading time4 min
Views4.4K

У меня возникла идея, как можно расширить синтаксис 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

Tags:
Hubs:
+6
Comments43

Articles