У меня возникла идея, как можно расширить синтаксис 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
Комментарии (43)
askv
31.05.2025 00:23result[>i] = e[i][>j][>k]a[>j]b[>k];Лёгким движением руки... 6 слагаемых превращаются в 27...
kinh Автор
31.05.2025 00:23Матрица e, вообще говоря, состоит из констант. Здесь надежда на то, что компилятор сообразит, что будет, если выражение умножить на 0, 1 или -1.
kovserg
31.05.2025 00:23А что мешает использовать анонимные функции для подобных операций?
struct eab_t { int i,j,k,v; }; template<class F> void eab(F f) { f({0,1,2, 1}); f({1,0,2,-1}); f({1,2,0, 1}); f({2,1,0,-1}); f({2,0,1, 1}); f({0,2,1,-1}); } ... double r[3],a[3],b[3]; eab([&](eab_t e){ r[e.i]=e.v*a[e.j]*b[e.k]; });
nin-jin
31.05.2025 00:23Интересно, но я бы разделил свёртку от итерирования, чтобы можно было делать и поэлементые операции, да и свёртки бывает нужны разные, а не только суммирование.
И как предлагается перегружать этот оператор?
kinh Автор
31.05.2025 00:23я бы разделил свёртку от итерирования
У меня была мысль для чистого итерирования использовать ещё два новых оператора:
|> проход от начала к концу
<| проход от конца к началу
но это была смутная идея, и я её не прорабатывал. Тем более, что при той записи, которую я предложил в статье, вполне можно разделить итерирование и свёртку для стандартных случаев: если декларация индекса не повторяется в одночлене - это итерирование, а если мы декларировали его повторно в том же одночлене - свёртка.И как предлагается перегружать этот оператор?
Сам по себе оператор не перегружается. В принципе, обычным образом перегружается обращение к элементу массива по его индексу (operator []), и перегружаются операции сложения и умножения. Теоретически, можно было бы добавить какую-то специальную операцию умножения, которая на входе принимает два множества объектов, и что-то с ними делает, но это слишком всё усложнит.
У меня была основная идея в том, чтобы простым образом задавать умножение нестандартных матриц. То есть, теперь можно, как в последнем примере статьи, рассматривать массив объектов, у которых вектором-строкой является одно из полей, как полноценную матрицу.
qw1
31.05.2025 00:23Вопрос-то не в этом был.
У вас в записи скалярного произведенияfloat result = a[>] * b[>];умножение явно прописано, а сложение подразумевается.
Это плохо.
Почему там именно сложение, а не какая-то другая операция, задаваемая пользователем?
Только ради того, что сейчас у вас именно такая задача?kinh Автор
31.05.2025 00:23Почему там именно сложение
Я стремился к лаконичности синтаксиса, а это - наиболее стандартный случай.
а не какая-то другая операция, задаваемая пользователем?
Ну, можно, например, так усложнить запись:
float result = a[>] *:+ b[>];Сначала мы поэлементно выполняем умножение, затем то, что получилось - складываем между собой.
nin-jin
31.05.2025 00:23Тогда не понятно как использовать произвольные функции вместо операторов. Ну а перегрузки нужны, чтобы иметь тот же лаконичный синтаксис и со своими типами данных.
kinh Автор
31.05.2025 00:23Тогда не понятно как использовать произвольные функции вместо операторов
Было бы проще, если бы вы привели пример, где это нужно.
Сформулируем задачу: сначала применяем поэлементно к двум векторам какую-то бинарную функцию, затем все элементы вектора результата обрабатываем другой функцией, аналогичной сложению. Для этого можно придумать такую запись:
class Obj {/*...*/}; Obj mul(Obj x, Obj y) {/*...*/} Obj add(Obj x, Obj y) {/*...*/} Obj a[10], b[10]; // ... Obj result = a[>] (mul:add) b[>];kinh Автор
31.05.2025 00:23Если нужно найти максимальное значение вектора, то это можно было бы так записать:
float a[10]; // ... float result = a[>] (: [](float x, float y) {return x > y? x: y;});Здесь внутри (: ...) мы задаём функцию от двух аргументов, которая сначала применяется к двум первым элементам вектора, затем её результат будет использоваться как левый аргумент этой же функции на следующих итерациях, а правый будет очередным элементом вектора. И так до конца вектора.
Если вектор имеет размер равный 1, то возвращается его первый элемент. Для учёта вектора нулевого размера можно дополнительно указать результат для этого случая после ключевого слова else.
Например, сумма всех элементов вектора:float a[10]; // ... float result = a[>] (: [](float x, float y) {return x + y;} else 0);qw1
31.05.2025 00:23Уже лучге, можно было бы записать скалярное произведение
float result = (a[>i] * b[>i]) (: [](float x, float y) {return x + y;} );Разделение операций по конвейеру (сначала делаем одно, потом результат передаём другому) очень похоже на ranges в новом стандарте.
Жаль конечно, что синтаксис лямбд в C++ не смогли сделать лаконичным.
Во многих языках уже ввели стрелочную нотацию(x,y) => x+yс авто-выводом типов параметров и результата. В таком виде ranges были бы очень хорошо читаемыми.nin-jin
31.05.2025 00:23В D можно вообще так писать:
list.map!q{ a + b }qw1
31.05.2025 00:23А почему конкретно a и b? Там такое соглашение, что параметры безымянной функции - это a, b, c, ... ? А если переменная b объявлена во внешней области видимости, есть способ к ней обратиться?
nin-jin
31.05.2025 00:23Тут фактически передаётся строка времени компиляции, так что ко внешним переменным обратиться нельзя. Чтобы было можно, нужно уже передавать полноценный делегат:
list.map!( a => a + b )
kinh Автор
31.05.2025 00:23В D можно вообще так писать: list.map!q{ a + b }
В моих обозначениях можно даже короче записать:
auto sum = list[>] (:+);Вместо бинарной функции мы передаём бинарный оператор.
nin-jin
31.05.2025 00:23На самом деле в D можно и так реализовать без изменения компилятора:
list.fold!q{+}
kinh Автор
31.05.2025 00:23Тогда и нет смысла как-то расширять тот синтаксис, что я описал в статье. В языке D остальные случаи гораздо более проработаны.
MaxAkaAltmer
31.05.2025 00:23Потом кто это месиво читать будет?
Создайте класс для подобных операций и все у вас будет хорошо.
KaminskyIlya
31.05.2025 00:23Даю голову на отсечение, что и сам автор лет через 30, глянув на свой же код будет долго рвать волосы на голове и других местах в попытках понять, что же все-таки в этом коде происходит. Это сейчас: да... по горячим следам, автору все кажется очевидным и само-собой разумеющимся (по себе сужу, да - я уже стар )
Вопрос к автору: чем вам не угодила перегрузка операторов? Вы можете для класса вектор перегрузить операцию - умножение на матрицу и писать примерно так:
Vector a = new Vector(1,2,3); Matrix m = new Matrix(1, 0, 1 2, 3, 2 1, -1, 4); Vector c = a*m; // умножение матриц Matrix scaling = new Matrix(...); Matrix translate = new Matrix(...); Matrix rotate = new Matrix(...); // деклаем матрицу, котора сразу и масштабирует и переносит и вращает Matrix transform = scaling * translate * rotate; // считаем определитель матрицы double det = transform;Так будет более понятно и не только вам, но и другим читателям кода. А все потому, что многие уже знают (со школы или университета), как работает умножение вектора на матрицу, и какой в итоге должен быть результат.
Нравится функциональный стиль? - пишите на соответствующем языке. На Haskel, например, можно так загнаться, что вы уже и через неделю не сможете разобраться в своем коде.
kinh Автор
31.05.2025 00:23Даю голову на отсечение, что и сам автор лет через 30, глянув на свой же код будет долго рвать волосы на голове и других местах в попытках понять, что же все-таки в этом коде происходит.
Запись полностью аналогична записи в тензорных обозначениях, где по повторяющимся индексам подразумевается суммирование. В физике такую запись используют как раз потому, что она удобнее и понятнее, чем безиндексная.
чем вам не угодила перегрузка операторов?
На Хабре уже была хорошая статья на тему того, почему умножение матриц без индексов - слишком неудобная операция: https://habr.com/ru/articles/910834/
Мне интересно, как вы будете умножать матрицы не строка на столбец, а строка на строку? А если строка в поле объекта, как у меня в последнем примере? Каждый раз будете оператор переопределять? Больше кода = больше ошибок.Нравится функциональный стиль?
В функциональном стиле там только переворачивание вектора. Я добавил это просто как интересную возможность использования прохода в обратном направлении. Это вполне можно выбросить - на умножение это не повлияет.
askv
31.05.2025 00:23Запись полностью аналогична записи в тензорных обозначениях, где по повторяющимся индексам подразумевается суммирование.
Только там она не называется скалярным произведением, это свёртка. Зачем Вы тогда её так называете?
Вообще для работы с матрицами есть специализированное ПО, тот же Матлаб.
kinh Автор
31.05.2025 00:23Только там она не называется скалярным произведением, это свёртка. Зачем Вы тогда её так называете?
Тут любое определение будет неточным.
Свёртка - это когда у одного и того же тензора два индекса одинаковы, и по ним производится суммирование. Скалярное произведение формально состоит из двух операций: сначала мы находим объект, состоящий из произведений каждого компонента одного вектора на каждый компонент второго, после чего проводим свёртку по обоим индексам этого объекта.Вообще для работы с матрицами есть специализированное ПО, тот же Матлаб.
Мне просто пришла интересная идея в голову, и я решил ей поделиться. В язык C++ это, разумеется, никто вносить не будет - но на Хабре хватает людей, которые пишут свои языки программирования, и им это может быть интересно.
askv
31.05.2025 00:23Скалярное произведение — это свёртка с метрическим тензором. Общая формула c=Gij u^i v^j. Только в ортонормированном базисе скалярное произведение выражается как сумма попарных произведений координат. Если базис другой, будет положительно определённая билинейная форма.
kinh Автор
31.05.2025 00:23Скалярное произведение — это свёртка с метрическим тензором
В моих обозначениях это так запишется:
auto c = g[>i][>j] * u[>i] * v[>j];Правда, отличать ковариантные и контравариантные компоненты тензора не получится.
KaminskyIlya
31.05.2025 00:23object IsHCons1 { type Aux[L[_], FH[_[_]], FT[_[_]], H0[_], T0[_] <: HList] = IsHCons1[L, FH, FT] { type H[t] = H0[t] ; type T[t] = T0[t] } def apply[L[_], FH[_[_]], FT[_[_]]](implicit tc: IsHCons1[L, FH, FT]): Aux[L, FH, FT, tc.H, tc.T] = tc implicit def mkIsHCons1[L[_], FH[_[_]], FT[_[_]]]: IsHCons1[L, FH, FT] = macro IsHCons1Macros.mkIsHCons1Impl[L, FH, FT] }Вот такой, например, понятный код. источник https://habr.com/ru/articles/259841/
sci_nov
31.05.2025 00:23Синтаксис Matlab содержит два умножения, типа а.*б и типа а*б, поэлементное и векторное. Может, проще сделать оператор точка-звездочка?
impwx
31.05.2025 00:23Кажется, что предлагаемый синтаксис подходит только для одной очень узкой задачи. Шаг в сторону - и, по всей видимости, нужно всё переписывать на обычные циклы. Например, если я хочу шаг, отличный от 1? А если я хочу на каждом шаге цикла еще залогировать что-то? Ну и так далее
goldexer
31.05.2025 00:23В каждом языке есть конструкции, о которых со временем понимают, что можно было организовать для них какой-то синтаксический сахар. А ведь когда-то давно в девяностых во всех разновидностях Basic-а или Pascal/Delphi очень не хватало «i++» инкремента и приходилось лепить конструкции вида «i=i+1», что прям сильно подбешивало. А ещё копировать массивы в цикле, тогда как весь цивилизованный мир использовал указатели/передачу по ссылке. Нет, ну способы были, конечно, но это небыло нативно, лаконично, удобно... Хорошо хоть строку - в виде массива символов - не надо было присваивать вручную по-одному... Если какие-то примитивы и инструкции уже оптимизированы и отработаны на уровне языка, то почему бы и нет? Оптимизация рутинных действий - итерировать, передать по цепочке, обработать базовые структуры - я за! Пока что такие операции вываливаются в изобретение велосипедов и эти самые велосипеды будут защищать, вместо того, чтоб от них избавляться.
HemulGM
31.05.2025 00:23В Delphi всегда был inc(I) и массивы всегда передавались по ссылке. А строки в Delphi всегда были куда проще в использовании, чем в С/Срр.
Nemoumbra
31.05.2025 00:23Если точнее, то inc появился ещё в турбо-паскале. У строк интрерфейс был удобный, но с нюансом про короткую длину строки.
qw1
31.05.2025 00:23Ещё точнее, это Вирт придумал.
Он загонялся по математической строгости, пытался различать ординалы и кардиналы. У него функции succ/prev и соответствующие им процедуры inc/dec, это не просто прибавить/убавить единицу, а выбрать предыдущий/следующий элемент в упорядоченном множестве, не обязательно числовом.
HemulGM
31.05.2025 00:23В те времена длины строки хватало, а сейчас длина до 2гб. И в редакторе кода можно в переменную вручную затолкать пару тысяч символов.
Apoheliy
31.05.2025 00:23По-моему, в текущем виде это всё очень притянуто за уши и слишком "местечково": если слева от "=" - то итерируем, если справа от "=" то что-то итерируем-складываем. Т.е. в одном случае одномерное остаётся одномерным, а в другом случае одномерное становится скаляром.
Возможно, выходом было бы отделить итерирование от всего остального. Например, ваш пример:
C[>i][>j] = A[i][>k] * B[>k][j];Превратить во что-то:
C[>i][>j] = SUMM(A[i][>k] * B[k][j]);Какие отличия:явное задание количества "итерирований": есть по i, j, k - по одной штуке. У
B[k][j] своей итерации нет;результат
A[i][>k] * B[k][j] - это одномерный массив (так как всего одна итерация по k) из результатов умножений;функция суммирования SUMM задаётся явно, принимает массив чего-то, возвращает сумму. Технически сюда можно и пользовательские функции затолкать.
Как результат, получаем в целом компактную запись с отработкой всяких случаев, когда элементами массивов являются разные штуки (условно: структуры, которые можно умножать и складывать).
Также, отдельно лучше/придётся продумать, когда индексы не-числовые: например, когда значениями индекса являются строки (или другие абстрактные ключи).
kovserg
31.05.2025 00:23Вся проблема в том что вы хотите использовать выражения для вычислений. А ведь можно просто разбить выражение на функции и обернуть их в функцию и её уже использовать по месту. Примерно так:
new: float result = a[>i] * b[>i]; old: float result = smul(a,b); new: a[>i] = i * i; old: for(auto i:index(a)) a[i]=i*i; new: a[>] = 0; old: load(a,0); new: result[>i] = a[i] * b[i]; old: mul(r,a,b); new: a[>i] = a[<i]; old: reverse(a);и не понадобиться изобретать инопланетный синтаксис к уже существующему.
kinh Автор
31.05.2025 00:23Я добавил в конец статьи ссылку на статью о проблемах безындексной записи умножения тензоров: Я не люблю NumPy
RainM
А в чем отличие от C/C++ Extensions for Array Notations? Это расширение языка было взято из Intel Cilk и поддерживалось когда-то классическим интеловским компилятором.
kinh Автор
Я не имел дела с этим. Сейчас посмотрел, и там нет суммирования по немым индексам, а пример умножения массивов выглядит так:
Насколько я понимаю - это не скалярное произведение, а поэлементное.