Краткий курс компьютерной графики: пишем упрощённый OpenGL своими руками, статья 3.14 из 6

  • Tutorial

Содержание основного курса




Общение вне хабра

Если у вас есть вопросы, и вы не хотите задавать их в комментариях, или просто не имеете возможности писать в комментарии, присоединяйтесь к jabber-конференции 3d@conference.sudouser.ru

4 Приветствие и вступление

Нумерация в прошлой статье закончилась на 3, в этой будем продолжать нумеровать насквозь.
UPD: ВНИМАНИЕ! Раздел, начиная с номера 3.1, 3.14 и 3.141 и далее, будет о тонкостях реализации основы основ компьютерной графики — линейной алгебры и вычислительной геометрии. О принципах графики пишет haqreu, я же буду писать о том, как это можно внятно запрограммировать!

Эта статья является продолжением серии статей о практической реализации элементов вычислительной геометрии, и, в частности, программного отрисовщика, с использованием C++98. Мы с haqreu сознательно идем на использование прошлой версии стандарта и написание собственной геометрической библиотеки для того, чтобы, во-первых, выпустить код примеров, которые без особых трудностей будут компилироваться большинством имеющихся компиляторов, а во-вторых, чтобы в нашем коде не было ничего, что скрыто в недрах библиотеки. В статье излагаются вопросы реализации шаблона прямоугольной матрицы template<size_t DimRows,size_t DimCols,typename number_t> class mat;

4.1 Благодарности
Я выражаю огромную признательность haqreu, как основоположнику данного курса. Так держать!
Я очень признателен lemelisk за предварительное рецензирование и ревью моих исходников. Спасибо за плодотворные дискуссии!
Также я должен поблагодарить Mingun за ценное замечание об оформлении шаблонов. Надеюсь, они стали доступнее для прочтения.

4.2 Небольшое философское замечание о языке, математике, C++, и всем таком
В целом, такие абстракции, как «матрица» и «вектор» и операции над ними необходимы для того, чтобы упростить и формализовать язык изложения геометрии. Без них нам приходилось бы постоянно говорить (и писать в программе) что-то вроде:

Говорим:
Чтобы получить новые координаты точки Point, которая перемещена в направлении Direction, нужно сложить
x точки Point с x перемещения Direction,
у точки Point с y перемещения Direction,
z точки Point с z перемещения Direction
Пишем:
Point.x=Point.x+Direction.x;
Point.y=Point.y+Direction.y;
Point.z=Point.x+Direction.z; //(Привет, "Эффект последней строки")

После того, как понятие «вектор» вошло в наш речевой оборот, мы уже можем говорить:
Чтобы получить новый радиус-вектор точки Point после перемещения в направлении Direction, нужно сложить векторы Point и Direction

А после включения в наш код шаблона vec, мы можем писать:
Point=Point+Direction;
Таким образом, мы сокращаем написанное, страхуемся от ошибок при копировании-вставке-правке и получаем способ рассуждать более емкими терминами.

Вот что об этом пишет Джефф Элджер в книге «C++ for Real programmers», (издательство AP Professional, перевод выпущен издательством «Питер» под названием «C++: Библиотека программиста»)

С++ —это на самом деле не столько язык, сколько инструмент для создания ваших собственных языков. Его элегантность заключается отнюдь не в простоте (слова С++ и простота режут слух своим явным противоречием), а в его потенциальных возможностях. За каждой уродливой проблемой прячется какая-нибудь умная идиома, изящный языковой финт, благодаря которому проблема тает прямо на глазах. Проблема решается так же элегантно, как это сделал бы настоящий язык типа Smalltalk или Lisp, но при этом ваш процессор не дымится от напряжения, а на Уолл-Стрит не растут акции производителей чипов памяти

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

5 Шаблон класса для обработки и хранения матриц

Мы храним матрицу как массив вектор-строк, используя соответствующую специализацию шаблона vec. Наш шаблон ориентирован только на работу с матрицами малой размерности (максимум 10x10 элементов. Хотя это уже будет ужас). Но в вычислительной геометрии, как правило, вся работа идет именно с матрицами не больше, чем 4x4. Большие же матрицы получаются, как правило, очень разреженными — для них заранее известно, что они практически полностью заполнены нулями, что позволяет оптимизировать их обработку и хранение.

5.1 Индексные операторы
Мы определили константный и неконстантный варианты индексного оператора [] таким образом, чтобы они возвращали ссылку на соответствующую вектор-строку матрицы:
их исходные тексты
Константный вариант
Неконстантный вариант
const vec<DimCols,number_t>& operator[] (const size_t idx) const 
{
        assert(idx<DimRows);
        return rows[idx];
}
vec<DimCols,number_t>& operator[] (const size_t idx) 
{
        assert(idx<DimRows);
        return rows[idx];
}
Два варианта операторов нужны из-за того, что мы везде, где только возможно, применяем const. Если бы не было константного варианта индексного оператора, компилятор не позволял бы нам собрать код, в котором, например, оператор вызывается по константной ссылке на экземпляр матрицы:
Так ругается GCC 4.9.2, когда не находит константного варианта оператора []
geometry.h:182: ошибка: no match for 'operator[]' (operand types are 'const mat<2ul, 3ul, float>' and 'size_t {aka long unsigned int}')
     for (size_t i=DimRows; i--; ret[i]=lhs[i]*rhs);
                                           ^
Кроме того, мы здесь применяем макрос assert() для отлова попыток обратиться к несуществующей строке матрицы. Подробнее про assert написано в статье 3.1

5.2 Вычисление обратной матрицы
Нам потребовалась подпрограмма для вычисления обратной матрицы. В частности, используя обратную матрицу можно «красиво» находить барицентрические координаты (можно посмотреть в моей ветке на гитхабе, функция filltria).
Задача вычисления обратной матрицы в общем случае является достаточно затратной по времени и может сопровождаться накоплением значительной погрешности при вычислениях. Однако для столь малых матриц эти проблемы не так существенны (если не строить такой пример специально). Мы выбрали алгоритм с применением союзной матрицы, потому как он показался нам более прозрачным, нежели алгоритм Гаусса-Жордана. Используется следующий факт:
  • если для данной матрицы составить союзную матрицу,
  • а затем транспонировать
  • и поделить ее на определитель данной,
мы получим обратную к матрицу. Это всего три внятных шага.
Союзная же матрица составляется из алгебраических дополнений. Алгебраическое дополнение — это определитель подматрицы, полученной из основной путем удаления строки и столбца, содержащих элемент, для которого мы находим это дополнение.
Пояснительная картинка из Википедии, автор - Александр Мехоношин
image

Получили следующий вопрос — как найти определитель.
Используем такое свойство:
определитель равен произведению элементов строки или столбца матрицы на их алгебраические дополнения.

Получилась рекурсия — алгебраическое дополнение — тоже определитель матрицы, только размерность на единицу меньше. А определитель матрицы размерности 1 — это ее единственный элемент. Рекурсия замкнулась. Мы воспользуемся прекрасной особенностью шаблонов C++ — развернем эту рекурсию прямо во время компиляции.

Замечание: Матрица у нас — прямоугольная, а определитель можно вычислить только для квадратной матрицы. Хотелось бы уже на этапе компиляции отлавливать попытки найти определитель неквадратной матрицы. Есть два способа это сделать:
  • унаследовать от mat класс squareMat, в котором объявить определитель
  • воспользоваться тем, что размеры матрицы нам известны еще на этапе компиляции, поэтому мы можем их сравнить и прервать компиляцию, если они не совпали

В первом случае мы получим отдельный тип матриц и кучу проблем с преобразованием матриц туда-сюда, потому как у нас может образоваться квадратная матрица в рамках типа mat, но чтобы посчитать определитель, нам придется конструировать из нее экземпляр типа squareMat.
Поэтому мы выберем второй способ — при разворачивании шаблонов обнаруживать несовпадение размеров и ломать компиляцию. В С++11 есть стандартное средство это сделать — static_assert(). В С++98 такого средства нет, поэтому мы изобретем суррогатное:
template<bool>
struct static_assertion;

template<>
struct static_assertion<true> {};
Если в этот шаблон подставлять true, просто образуется пустая структура static_assertion{}. А вот если поставить false, образуется ошибка (попытка использовать неполное объявление структуры static_assertion), что и останавливает компиляцию:
geometry.h:125: ошибка: invalid use of incomplete type 'struct static_assertion<false>'
          static_assertion<DimCols<=10>();
          ^
Чем больше ошибок смещено на этап компиляции, тем лучше, ведь они просто не допускают плохой код до работы.
Вернемся к нашей рекурсии:
Исходные тексты шаблона, отвечающего за рекурсию
/////////////////////////////////////////////  полный шаблон структуры dt
1.   template<size_t Dim,typename number_t> struct dt 
2.   {
3.      static number_t det(const mat<Dim,Dim,number_t>& src) 
4.      {
5.         number_t ret=0;
6.         for (size_t i=Dim; i--;)
7.         {
8.               ret += src[0][i] * src.cofactor(0,i);
9.         }
10.       return ret;
11.    }
12. };
/////////////////////////////////////////////  частичная специализация шаблона для случая, когда размерность равна 1.
13. template<typename number_t> struct dt<1,number_t> 
14. {
15.     static number_t det(const mat<1,1,number_t>& src) 
16.     {
17.         return src[0][0];
18.     }
19. };


Итак, для всех ситуаций, когда первый параметр шаблона не равен 1, будет применяться верхний шаблон (строки 1-11). Для ситуации, первый параметр равен 1, определена специальная версия шаблона (строки 12-18). Обратите внимание на строку 13: Дополнение <1,number_t> как раз и указывает, что это объявление — отдельная сущность. С общей версией шаблона dt ее ничего не связывает. Встречается ошибка, когда путают поведение частичной специализации и наследования, ожидая, что все объявленное в общей версии окажется и в частичной. Этого не будет — на программиста возлагается воссоздание в частичной специализации всего, что ему нужно.
  • Оборачивание в отдельную структуру выполнено потому, что стандарт не разрешает делать частичную специализацию для функций.
  • Мы вынесли эту функцию из mat потому, что поступив иначе, были бы вынуждены дописывать в mat<1,1,number_t> часть методов из общей версии, что привело бы к некоторому раздутию кода.
Частичная специализация, точно по определению определителя матицы 1x1, отдает нам единственный элемент (строка 16).

Общая же версия шаблона выполняет разложение матрицы по нулевой строке — она берет элемент с индексом 0, i из первой строки, умножает его на алгебраическое дополнение — cofactor(0,i) а результат накапливает в переменной ret. Все в полном согласии с формулой этого разложения:image

5.3 Как происходит рекурсия
ВАЖНО: Строго говоря, механизм обработки шаблонов — личное дело компилятора. Ниже я рассказываю о том, как это могло бы происходить. Это делается только для иллюстрации того, как рекурсивно разворачиваются декларации шаблонов.

Представим, что мы самостоятельно выполняем разворачивание шаблонов. Написан код:
...
using namespace std;
int main()
{
    mat<3,3,int> m;
    ....
    cout << m.det() << endl;
    return 0;
}
Мы встречаем в нем вызов det() из шаблона mat, а значит, должны получить текст функции det() для конкретного случая mat<3,3,int>.
Шаблон
Результат подстановки DimCols=3; DimRows=3; number_t=int
number_t det() const 
{
        static_assertion<DimCols==DimRows>();
        return dt<DimCols,number_t>::det(*this);
}
int det() const<3,int> 
{
        static_assertion<3==3>();
        return dt<3,int>::det(*this);
}
Отлично, теперь нам нужно из шаблона template<size_t Dim,typename number_t> struct dt получить конкретный текст для случая dt<3,int>:
Шаблон
Результат подстановки Dim=3; number_t=int
template<size_t Dim,typename number_t> struct dt 
{
   static number_t det(const mat<Dim,Dim,number_t>& src) 
   {
        number_t ret=0;
        for (size_t i=Dim; i--;)
        {
              ret += src[0][i] * src.cofactor(0,i);
        }
       return ret;
    }
};
struct dt<3,int>
{
     static int det(const mat<3,3,int>& src) 
     {
         int ret=0;
         for (size_t i=3; i--;)
         {
               ret += src[0][i] * src.cofactor(0,i);
         }
       return ret;
    }
};
Мы встретили вызов оператора [] и cofactor() из mat<3,3,int>. Тело [] сейчас для нас интереса не представляет, разбираемся с cofactor():
Шаблон
Результат подстановки DimCols=3; DimRows=3; number_t=int
number_t cofactor(size_t row, size_t col) const 
{
        return get_minor(row,col).det()*((row+col)%2 ? -1 : 1);
}
int cofactor(size_t row, size_t col) const 
{
        return get_minor(row,col).det() *( (row + col) % 2 ? -1 : 1);
}
Здесь вызывается get_minor(). Получим для нее текст:
Шаблон
Результат подстановки DimCols=3; DimRows=3; number_t=int
mat<DimRows-1,DimCols-1,number_t> get_minor(size_t row, size_t col) const 
{
     mat<DimRows-1,DimCols-1,number_t> ret;
     for (size_t i=DimRows-1; i--; )
     {
            for (size_t j=DimCols-1;j--;)
            {
                   ret[i][j]=rows[ i<row ? i : i+1 ] [  j <col ? j : j+1 ];   
// конструкция j <col ? j : j+1 обеспечивает пропуск столбца с номером col
//аналогично  i<row ? i : i+1 - обеспечивает пропуск строки с номером row
            }                                                                     
     }
     return ret;
}
mat<2,2,int> get_minor(size_t row, size_t col) const 
{
     mat<2,2,int> ret;
     for (size_t i=2; i--; )
     {
            for (size_t j=2;j--;)
            {
                   ret[i][j]=rows[ i<row?i:i+1 ] [ j<col?j:j+1 ];
            }
     }
     return ret;
}
Мы получили, что этот get_minor() возвращает уже матрицу 2x2. Теперь мы, для матрицы 2x2 будем вызывать mat<2,2,int>::det(). Но у нас есть только тело mat<3,3,int>::det(). Поэтому мы совершаем погружение в рекурсию на один шаг и строим еще один набор методов:
Методы, появившиеся при разворачивании шаблона 2x2
Мы должны получить текст функции det() для конкретного случая mat<2,2,int>:
Шаблон
Результат подстановки Dim=2; number_t=int
number_t det() const 
{
        static_assertion<DimCols==DimRows>();
        return dt<DimCols,number_t>::det(*this);
}
int det() const<2,int> 
{
        static_assertion<2==2>();
        return dt<2,int>::det(*this);
}
Отлично, теперь нам нужно из шаблона template<size_t Dim,typename number_t> struct dt получить конкретный текст для случая dt<2,int>:
Шаблон
Результат подстановки Dim=2; number_t=int
template<size_t Dim,typename number_t> struct dt 
{
   static number_t det(const mat<Dim,Dim,number_t>& src) 
   {
        number_t ret=0;
        for (size_t i=Dim; i--;)
        {
              ret += src[0][i] * src.cofactor(0,i);
        }
       return ret;
    }
};
struct dt<2,int>
{
     static int det(const mat<2,2,int>& src) 
     {
         int ret=0;
         for (size_t i=2; i--;)
         {
               ret += src[0][i] * src.cofactor(0,i);
         }
       return ret;
    }
};
Мы встретили вызов оператора [] и cofactor() из mat<2,2,int>. Тело [] сейчас для нас интереса не представляет, разбираемся с cofactor():
Шаблон
Результат подстановки DimCols=2; DimRows=2; number_t=int
number_t cofactor(size_t row, size_t col) const 
{
        return get_minor(row,col).det()*((row+col)%2 ? -1 : 1);
}
int cofactor(size_t row, size_t col) const 
{
        return get_minor(row,col).det() *( (row + col) % 2 ? -1 : 1);
}
Здесь вызывается get_minor(). Получим для нее текст:
Шаблон
Результат подстановки DimCols=2; DimRows=2; number_t=int
mat<DimRows-1,DimCols-1,number_t> get_minor(size_t row, size_t col) const 
{
     mat<DimRows-1,DimCols-1,number_t> ret;
     for (size_t i=DimRows-1; i--; )
     {
            for (size_t j=DimCols-1;j--;)
            {
                   ret[i][j]=rows[ i<row?i:i+1 ] [ j<col?j:j+1 ];
            }
     }
     return ret;
}
mat<1,1,int> get_minor(size_t row, size_t col) const 
{
     mat<1,1,int> ret;
     for (size_t i=1; i--; )
     {
            for (size_t j=1;j--;)
            {
                   ret[i][j]=rows[ i<row?i:i+1 ] [ j<col?j:j+1 ];
            }
     }
     return ret;
}


Развернув цепочку шаблонов, мы снова столкнулись с вызовом get_minor() в теле cofactor(), но уже для случая DimCols=2, DimRows=2. Полученная версия get_minor() возвращает матрицу 1x1. Далее, для нее в теле cofactor() запрашивается определитель. Мы снова должны получить текст mat::det(), но уже для случая DimCols=1; DimRows=1.
Шаблон
Результат подстановки DimCols=1; DimRows=1; number_t=int
number_t det() const 
{
        static_assertion<DimCols==DimRows>();
        return dt<DimCols,number_t>::det(*this);
}
int det() const<1,int> 
{
        static_assertion<1==1>();
        return dt<1,int>::det(*this);
}
Отлично, теперь нам нужно из шаблона template<size_t Dim,typename number_t> struct dt получить конкретный текст для случая dt<1,int>. СТОП, это же попадает под частичную специализацию шаблона dt<1,number_t>! Нужно обрабатывать именно ее, а не общую версию шаблона!
Шаблон
Результат подстановки Dim=1; number_t=int
template<typename number_t> struct dt<1,number_t> 
{
     static int det(const mat<1,1,number_t>& src) 
     {
         return src[0][0];
     }
};
template<typename number_t> struct dt<1,int> 
{
     static number_t det(const mat<1,1,int>& src) 
     {
         return src[0][0];
     }
};

Видим, что здесь уже нет обращения к cofactor(). На этом раскрутка рекурсии закончена — мы получили тела функций для вычисления определителя матриц 3x3, 2x2, 1x1.
Важное замечание: Может возникнуть желание написать так:
number_t det()
{
    if(DimCols==1 && DimRows==1)
    {
         return rows[0][0];
    }
    number_t ret=0;
    for (size_t i=DimCols; i--;)
    {
               ret += src[0][i] * src.cofactor(0,i);
    }
    return ret;
}

Мотивация простая — это будет вроде как обычная рекурсия, да и частичная специализация не нужна. Однако во время развертывания шаблонов, компилятор ничего не будет делать с условным оператором, а начнет строить шаблоны(3, 2, 1, 0, std::numeric_limits<size_t>::max, std::numeric_limits<size_t>::max-1, ...), пока его не остановит внутренний счетчик -ftemplate-depth=(по умолчанию, эта глубина погружения равна 900).
Возможность использовать привычные операторы управления потоком появляется в C++11 для функций, объявленных c ключевым словом constexpr. Они позволяют организовать еще один вид вычислений во время компиляции, отличный от рекурсии на шаблонах. А не провернуть ли мне весь наш рендер во время компиляции?

Итак, мы описали в терминах C++ отыскание определителя матрицы, а также сопутствующие этому операции — построение минора матрицы и алгебраического дополнения. Для отыскания обратной матрицы осталось только построить союзную матрицу:
Метод, строящий союзную матрицу:
mat<DimRows,DimCols,number_t> adjugate() const 
{
        mat<DimRows,DimCols,number_t> ret;
        for (size_t i=DimRows; i--;)
        {
            for (size_t j=DimCols; j--;)
            {
                 ret[i][j]=cofactor(i,j)
            }
        }    
        return ret;
}
действует точно по определению — заполняет матрицу алгебраическими дополнениями соответствующих элементов.
Осталось выполнить главный этап построений — найти обратную матрицу (она будет получена в транспонированном виде):
mat<DimRows,DimCols,number_t> invert_transpose() 
{
        mat<DimRows,DimCols,number_t> ret = adjugate();  //получили союзную матрицу

        number_t det = ret[0]*rows[0];                   //вычислили определитель исходной матрицы, путем скалярного
                                                  //умножения первой вектор-строки
                                                  //элементов на вектор строку их алгебраических дополнений

        return ret/det;                           //поделили союзную матрицу на определитель исходной,
                                                  //получив тем самым транспонированную обратную матрицу
}

5.4 Умножение матриц
Операция перемножения матриц также довольно затратна по времени и может сопровождаться накоплением погрешности. В случаях матриц большей размерности (больше, чем 64x64), как правило, применяют специальные алгоритмы, вроде алгоритма Штрассена. Асимптотический рекордсмен — Алгоритм Копперсмита-Винограда (с улучшениями Вирджинии Вильямс), сейчас не применяется, потому как начинает обгонять алгоритм Штрассена на матрицах настолько огромных, что в современной вычислительной практике еще не встречаются. У нас матрицы мелкие, поэтому нам для умножения подходит просто определение операции умножения матриц:
Произведение матриц AB состоит из всех возможных комбинаций скалярных произведений вектор-строк матрицы A и вектор-столбцов матрицы B

image
Сведения об авторе картинки
«Matrix multiplication diagram 2» участника File:Matrix multiplication diagram.svg:User:BilouСм. ниже — Этот файл является производной работой от: Matrix multiplication diagram.svg. Под лицензией CC BY-SA 3.0 с сайта Викисклада.

Определение в таком виде позволяет нам использовать операцию скалярного умножения из нашего шаблона vec. Но для этого нам нужен метод, извлекающий из матрицы ее столбец в виде вектора:
текст метода col()
vec<DimRows,number_t> col(const size_t idx) const 
{
        assert(idx<DimCols);
        vec<DimRows,number_t> ret;
        for (size_t i=DimRows; i--;)
        {
               ret[i]=rows[i][idx]; 
        }
        return ret;
}
Да, здесь вроде как имеет место быть лишнее копирование (позже мы посмотрим, справится ли компилятор с его выбрасыванием), но главная цель — написать такое краткое и емкое умножение матриц:
template<size_t R1,size_t C1,size_t C2,typename number_t>
    mat<R1,C2,number_t> operator*(
                                     const mat<R1,C1,number_t>& lhs, 
                                     const mat<C1,C2,number_t>& rhs
                                 ) 
{
    mat<R1,C2,number_t> result;
    for (size_t i=R1; i--; )
    {
        for (size_t j=C2; j--;)
        {
              result[i][j]=lhs[i]*rhs.col(j);               //перемножить скалярно i-тую строку левого операнда
                                                            // на j-тый столбец правого
        }
    }
    return result;
}

Уделим внимание этому тексту. Здесь на C++ изложено математическое определение операции умножения прямоугольных матриц, причем сделано это очень близко к тексту[«заполнение матрицы-результата всевозможными скалярными произведениями вектор-строк левой матрицы на вектор-столбцы правой»].
При этом контроль размерностей (число столбцов левой должно быть равно числу строк правой) выполняется сразу на этапе компиляции — попытка перемножить не подходящие по размеру матрицы даже не откомпилируется.
А это очень важно, ведь ошибка времени выполнения может быть плавающей, трудноуловимой и так далее. А ошибка компиляции вылезает практически всегда. (Нужно учесть сумасшествие IDE, которая может часть проекта закэшировать, часть забыть обновить и так далее. Но от этого есть верное лекарство — Rebuild All).

5.5 Генератор единичных матриц
В его тексте:
static mat<DimRows,DimCols,number_t> identity() 
{
        mat<DimRows,DimCols,number_t> ret;
        for (size_t i=DimRows; i--; )
        {
            for (size_t j=DimCols;j--;)
            {
                   ret[i][j]=(i==j) ;
            }
        }
        return ret;
}

Есть интересная особенность — запись результата типа bool в переменную априорно неизвестного типа number_t. Здесь может случиться кошмарик у человека, который захотел бы использовать нашу геометрию в своем проекте, причем вместе с самодельным типом для хранения чисел (это может быть тип, который хранит числа с фиксированной запятой). Он может не ожидать, что наши методы будут пытаться записывать в его тип значения bool (оператор == для size_t возвращает bool), и нагородить собственный конструктор, принимающий bool и делающий это преобразование не так, как описано в стандарте ([§4.7/4] требует, чтобы при приведении bool в числовой тип, false превращалось в 0, true — в 1). Так как в данном случае, от стандарта отступает гипотетический пользователь, нам не следует пытаться подстилать ему соломку (мы могли бы сотворить такое: ret[i][j]=static_cast<number_t>(static_cast<int>(i==j)), и передавать ему хотя бы int) ценой череды преобразований — нарушил стандарт — сам виноват.

5.6 Умножение матрицы на вектор
Для умножения матрицы (слева) на вектор (справа), мы рассмотрим вектор как матрицу из одного столбца. Далее, мы можем выполнить уже известную нам операцию умножения матриц, которая в итоге также даст матрицу из одного столбца — этот вектор-столбец и будет результатом умножения.
Нам понадобится метод set_col() для заполнения столбца матрицы из вектора:
void set_col(size_t idx, vec<DimRows,number_t> v) 
{
        assert(idx<DimCols);
        for (size_t i=DimRows; i--;)
        {
                rows[i][idx]=v[i];
        }
}

И генератор матрицы по заданному вектор-столбцу
static mat<DimRows,1,number_t> make_from_vec_col(const vec<DimRows,number_t>& v)
{
        static_assertion<DimCols==1>();
        mat<DimRows,1,number_t> ret;
        ret.set_col(0,v);
        return ret;
}

Теперь мы можем записать оператор умножения матрицы на вектор:
template<size_t DimRows,size_t DimCols,typename number_t> 
    vec<DimRows,number_t> operator*(
                                        const mat<DimRows,DimCols,number_t>& lhs, 
                                        const vec<DimCols,number_t>& rhs
                                   )
{
    return (lhs * mat<DimCols,1,number_t>::make_from_vec_col(rhs) ).col(0);
}
Сработано точно по определению — получили из вектора матрицу с одним столбцом, перемножили, вернули столбец-результат.

5.7 Деление матрицы на скаляр, вывод в ostream
Выполняются с использованием того факта, что матрица у нас хранится как массив вектор — строк, а для вектора аналогичные операции уже написаны:
Тексты операций
template<size_t DimRows,size_t DimCols,typename number_t>
   mat<DimCols,DimRows,number_t> operator/(
                                             mat<DimRows,DimCols,number_t> lhs, 
                                             const number_t& rhs
                                          ) 
{
    for (size_t i=DimRows; i--; )
    {
         lhs[i]=lhs[i]/rhs;
    }
    return lhs;
}

template <size_t DimRows,size_t DimCols,class number_t> 
   std::ostream& operator<<(   std::ostream& out, 
                               const mat<DimRows,DimCols,number_t>& m
                           ) 
{
    for (size_t i=0; i<DimRows; i++) 
    {
              out << m[i] << std::endl;
    }
    return out;
}


6 Заключение

В этой статье я постарался максимально доступно описать построение шаблона класса mat. В следующей статье я затрону вопросы оптимизации и ускорения данного кода на архитектурах x86_64 и ARM (на примере STM32F3). До скорой статьи!
Share post
AdBlock has stolen the banner, but banners are not teeth — they will be back

More
Ads

Comments 64

    +1
    Прежде чем выйдет следующая статья — однопоточный растеризатор, основанный на half-space (в этом коммите), типично раза в 4 раза медленнее растеризатора на сканлайнах — за счёт оверхеда по bounding box и трех счетчиков. Чтобы извлечь приличную выгоду из half-space, нужно разбивать экран на блоки и распараллеливать отрисовку по блокам, как это делается в GPU. Выглядит он, конечно, просто, но на практике такой простой код никуда не годится.
      +4
      Секунду, <<на практике>> нас вообще не интересует, это вынесено в место, про которое отдельно сказано, что для нас это чёрный ящик. Чем проще код, тем лучше. Разумеется, в реальной жизни это всё запрограммировано даже не прошивкой, а тупо кремнием. Но это не суть данного курса. Данный курс нацелен на то, чтобы показать, что все вычисления за красивыми картинками тривиальны и доступны любому (заинтересованному) школьнику.
        +1
        >> Чтобы извлечь приличную выгоду из half-space, нужно разбивать экран на блоки и распараллеливать отрисовку по блокам, как это делается в GPU.

        Интересная статья про софтверную растеризацию в Larrabee.
        software.intel.com/en-us/articles/rasterization-on-larrabee
        0
        Хорошая статья, но есть одно но. Заголовки/текст и т.д. разных размеров не понятно зачем уж очень сильно режет глаз, особенно в начале статьи заметно. Хотя это лично мое имхо.
          +1
          Я использовал стандартный контейнер <h4> для заголовков верхнего уровня (4, 5, 6), и <h5> для заголовков вложенного (5.1, 5.2, ..., 5.7) уровня. Перед публикацией я посмотрел статью на «обычном» Firefox и Chrome, а также на Firefox с темой HabraDarkAges и не заметил, что текст визуально «рассыпается». В последнем случае, заголовки даже были выделены отдельным цветом, что на мой взгляд, удобно.
          +3
          Пожалейте же читателей кода! Зачем лепить километровые предложения template<...> на одну строку с самим типом/функцией? Ну проще же читать, когда можно разделить эти сущности визуально.
            0
            Я внес поправки и надеюсь, теперь шаблоны действительно стали выглядеть лучше.
            0
            Есть два способа это сделать:
            Справедливости ради, есть ещё третий способ: специализировать шаблон mat для квадратных матриц и добавить в эту специализация характерные только для таких матриц методы.
              0
              Здесь снова проблема: в эту частичную специализацию придется также поместить индексные операторы, хранение самой матрицы, генератор единичной матрицы и так далее. Удвоим объем кода: получим два места, в которых фактически делается одно и тоже.
                0
                Можно написать базовый класс с общими операциям и унаследовать оба варианта от него.
                  0
                  При наследовании шаблона от шаблона появятся дикие конструкции вроде this->rows[i][j]. Обсуждение этого феномена на stackoverflow
              0
              Я по поводу нахождения обратной матрицы. Можно обращать, используя теорему Гамильтона-Кэли. В этом случае нужно найти всего лишь n степеней матрицы и следы полученных матриц. Гораздо быстрее будет работать для больших матриц, чем вычисления с использованием союзной матрицы. Есть ещё метод Крылова.
                0
                Но ведь в самом начале текста сказано, что нас интересуют только маленькие плотнозаполненные матрицы…
                  0
                  Думаю, что стоит попробовать сравнить алгоритмы на матрицах малого порядка.
                    0
                    Было бы интересно, можете написать код?
                      +1
                      Попробую, только не тестировал (возможно, косякнул, прошу извинить). Сейчас кину ссылку на книгу, где описан метод. Вот код.

                      mat<DimRows,DimCols,number_t> invert_transpose(const mat<DimRows,DimCols,number_t>& A)
                      {
                      	mat<DimRows,DimCols,number_t> PowMat[DimRows+1];
                      
                      	// Вычисляем степени матриц
                      	PowMat[0] = identity();
                      	PowMat[1] = A;
                      	for (size_t i=2; i <= DimRows; i++)
                      	{
                      		PowMat[i] = PowMat[i-1] * A;
                      	}
                      
                      	number_t s[DimRows], p[DimRows];
                      	for (size_t i = 0; i < DimRows; i++)
                      	{
                      		// Вычисляем следы матриц (суммы диагональных элементов)
                      		s[i] = 0;
                      		for(size_t j = DimRows; j--; )
                      			s[i] += PowMat[i+1][j][j];
                      
                      		// Определяем коэффициенты характеристического уравнения матрицы A
                      		p[i] = s[i];
                      		for(size_t j = 1; j < i; j++)
                      			p[i] -= p[j-1] * s[i-j];
                      		p[i] /= i + 1;
                      	}
                      
                      	// Рассчитываем обратную матрицу по следствию из теоремы Гамильтона-Кэли
                      	mat<DimRows,DimCols,number_t> result = PowMat[DimRows - 1];
                      	for(size_t i = 0; i < DimRows - 1; i++)
                      		result -= p[i] * PowMat[DimRows - i - 2];
                      	return result/p[DimRows - 1];
                      }
                      
                        0
                        Здесь тип number_t предполагает либо double, либо long double, либо класс с перегруженными арифметическими операциями (например, mpreal).
                          0
                          Код работать не желает.
                          Исходная матрица:
                                                            [ 1  4  9 ]
                                                            [         ]
                          (%o4)                             [ 2  9  7 ]
                                                            [         ]
                                                            [ 4  3  8 ]
                          
                          

                          Решение при помощи maxima:
                          (%i3) float(transpose(invert(matrix([1,4,9],[2,9,7],[4,3,8]))));
                                [ - 0.2982456140350877  - 0.07017543859649122    0.1754385964912281   ]
                                [                                                                     ]
                          (%o3) [ 0.02923976608187134    0.1637426900584795    - 0.07602339181286549  ]
                                [                                                                     ]
                                [  0.3099415204678362   - 0.06432748538011696  - 0.005847953216374269 ]
                          

                          Решение нашим методом с союзной матрицей:
                          -0.298246 -0.0701754 0.175439 
                          0.0292398 0.163743 -0.0760234 
                          0.309942 -0.0643275 -0.00584795 
                          

                          Решение вашим методом:
                          0.352381 0.015873 0.168254 
                          -0.0380952 0.603175 -0.0349206 
                          0.0952381 -0.0412698 0.511111 
                          

                          Потребовалось доопределить отсутствующие у нас операторы - умножение матрицы на скаляр и вычетание матриц:
                          template<size_t DimRows,size_t DimCols,typename T> 
                              mat<DimRows,DimRows,T> operator-(
                                                             const mat<DimRows,DimCols,T>& lhs, 
                                                             const mat<DimRows,DimCols,T>& rhs
                                                 )
                          {
                              mat<DimRows,DimCols,T> ret;
                              for(size_t i=DimCols;i--;)
                              {
                                  ret[i]=lhs[i]-rhs[i];
                              }
                              return(ret);
                          }
                          template<size_t DimRows,size_t DimCols,typename T> mat<DimCols,DimRows,T>
                                  operator*(
                                                  mat<DimRows,DimCols,T> lhs, 
                                                  const T& rhs
                                                ) 
                          {
                              for (size_t i=DimRows; i--;)
                              {
                                  lhs[i]=lhs[i]*rhs;
                              }
                              return lhs;
                          }
                          
                            0
                            Я ещё раз сейчас проверю. Книга Воеводиных «Параллельные вычисления», стр. 204-205 + книга Гатмахера.
                              +1
                              У меня есть подозрение, что если просто сосчитать количество арифметических операций, потребных для применения этого метода, и сравнить с количеством операций, необходимых для метода с союзной матрицей — разница будет в пользу метода с союзной матрицей.

                              Кроме того, в этом учебнике бойко сообщают оценки числа потребных процессоров, но ничего не говорят об оценках потребности на межпроцессную коммуникацию: пересылку данных и синхронизацию. Это довольно странно.
                              0
                              *Гантмахера.
                              +2
                              Нашёл первую ошибку.

                              mat<DimRows,DimCols,number_t> invert_transpose(const mat<DimRows,DimCols,number_t>& A)
                              {
                                  mat<DimRows,DimCols,number_t> PowMat[DimRows+1];
                              
                                  // Вычисляем степени матрицы
                                  PowMat[0] = identity();
                                  PowMat[1] = A;
                                  for (size_t i=2; i <= DimRows; i++)
                                  {
                                      PowMat[i] = PowMat[i-1] * A;
                                  }
                              
                                  number_t s[DimRows], p[DimRows];
                                  for (size_t i = 0; i < DimRows; i++)
                                  {
                                      // Вычисляем следы матриц (суммы диагональных элементов)
                                      s[i] = 0;
                                      for(size_t j = DimRows; j--; )
                                          s[i] += PowMat[i+1][j][j];
                              
                                      // Определяем коэффициенты характеристического уравнения матрицы A
                                      p[i] = s[i];
                                      for(size_t j = 1; j <= i; j++)
                                          p[i] -= p[j-1] * s[i-j];
                                      p[i] /= i + 1;
                                  }
                              
                                  // Рассчитываем обратную матрицу по следствию из теоремы Гамильтона-Кэли
                                  mat<DimRows,DimCols,number_t> result = PowMat[DimRows - 1];
                                  for(size_t i = 0; i < DimRows - 1; i++)
                                      result -= p[i] * PowMat[DimRows - i - 2];
                                  return result/p[DimRows - 1];
                              }
                              

                              Учитывается равенство во внутреннем цикле. Попробуйте.
                                +1
                                Сошлось, поздравляю (учтите, что метод с союзной матрицей выдает транспонированное решение)!
                                ////////С союзной матрицей:
                                -0.298246 -0.0701754 0.175439 
                                0.0292398 0.163743 -0.0760234 
                                0.309942 -0.0643275 -0.00584795 
                                
                                ////////По теореме Кэли-Гамильтона
                                -0.298246 0.0292398 0.309942 
                                -0.0701754 0.163743 -0.0643275 
                                0.175439 -0.0760234 -0.00584795 
                                

                                Кстати тот факт, что метод с союзной матрицей «бесплатно» еще и транспонирует оказался нам на руку — с том месте, где нам нужно обращать матрицу, нам как раз нужен транспонированный вариант.
                                  0
                                  Не совсем, некоторые элементы получаются с ошибкой. А если потестировать ещё?
                                    +2
                                    Все одинаково, просто одну из матриц нужно транспонировать для сравнения
                                      0
                                      Всё понял, спасибо.
                                      +3
                                      Я вроде каждое с каждым сравнил (числа совпадают), сильно путает «транспонированность».

                                      Я сейчас задал сильно перекошенную матрицу (и транспонировал результат метода с союзной матрицей, чтобы глаза не ломать)
                                      Это вариант с float
                                      Source
                                                2.00000000           4.00000000           9.00000000           2.00000000 
                                                2.00000000           9.00000000           7.00000000      100000.00000000 
                                                4.00000000           3.00000000           8.00000000           0.00000100 
                                             4000.00000000        6633.00000000      995522.00000000           1.13000000 
                                      
                                      With Cofactors
                                               -0.29927686           0.00000599           0.40014544          -0.00000051 
                                                0.40298930          -0.00000806          -0.19945022          -0.00000204 
                                               -0.00148256           0.00000003          -0.00027888           0.00000102 
                                               -0.00003018           0.00001000           0.00000997           0.00000000 
                                      
                                      With Cayley–Hamilton theorem
                                               -0.30498824           0.00000610           0.40778178          -0.00000052 
                                                0.41067991          -0.00000824          -0.20325670          -0.00000872 
                                               -0.00151085           0.00000003          -0.00028420           0.00000104 
                                               -0.00003080           0.00001036           0.00000000           0.00000000 
                                      


                                      Это вариант с double
                                      Source
                                                2.00000000           4.00000000           9.00000000           2.00000000 
                                                2.00000000           9.00000000           7.00000000      100000.00000000 
                                                4.00000000           3.00000000           8.00000000           0.00000100 
                                             4000.00000000        6633.00000000      995522.00000000           1.13000000 
                                      
                                      With Cofactors
                                               -0.29927683           0.00000599           0.40014542          -0.00000051 
                                                0.40298926          -0.00000806          -0.19945022          -0.00000204 
                                               -0.00148256           0.00000003          -0.00027888           0.00000102 
                                               -0.00003018           0.00001000           0.00000997           0.00000000 
                                      
                                      With Cayley–Hamilton theorem
                                               -0.29927683           0.00000599           0.40014542          -0.00000051 
                                                0.40298926          -0.00000806          -0.19945022          -0.00000204 
                                               -0.00148256           0.00000003          -0.00027888           0.00000102 
                                               -0.00003018           0.00001000           0.00000997           0.00000000
                                      


                                      Это насчитала maxima
                                           [ - 2.992768282b-1   5.985538325b-6    4.001454189b-1   - 5.099975655b-7 ]
                                           [                                                                        ]
                                           [  4.029892624b-1   - 8.059760198b-6  - 1.994502193b-1  - 2.04038202b-6  ]
                                           [                                                                        ]
                                           [ - 1.482559332b-3   2.963966182b-8   - 2.788771941b-4   1.02014204b-6   ]
                                           [                                                                        ]
                                           [ - 3.01797179b-5    1.000060359b-5    9.967132764b-6   1.224243903b-10  ]
                                      

                                      Можно заметить, что на float вариант по теореме немного отстает. Но это в данном конкретном случае.
                                        +1
                                        А что со временем работы?
                                          +2
                                          Одна и та же матрица миллион раз обращалась.
                                             for(size_t i=0;i<1000000;i++)
                                              {
                                                  alpha=alpha.invert_transpose();
                                              }
                                          
                                          Вариант с союзной матрицей: 390 308 мкс
                                                    2.00000000           4.00000000           9.00000000           2.00000000 
                                                    2.00000000           9.00000000           7.00000000      100000.00000000 
                                                    4.00000000           3.00000000           8.00000000           0.00000100 
                                                 4000.00000000        6633.00000000      995522.00000000           1.13000000 
                                          390 308 микросекунд спустя ----------------------------
                                                    2.00000000           4.00000000           9.00000000           2.00000000 
                                                    2.00000000           9.00000000           7.00000000       99999.99999999 
                                                    4.00000000           3.00000000           8.00000000           0.00000100 
                                                 4000.00000000        6633.00000000      995521.99999991           1.13000000
                                          


                                          Вариант по теореме: 1 291 526 мкс
                                                     2.00000000           4.00000000           9.00000000           2.00000000 
                                                    2.00000000           9.00000000           7.00000000      100000.00000000 
                                                    4.00000000           3.00000000           8.00000000           0.00000100 
                                                 4000.00000000        6633.00000000      995522.00000000           1.13000000 
                                          1 291 526 микросекунд спустя----------------------
                                                    1.99916579           3.99984523           8.99719176           1.99992263 
                                                    2.00265273           9.00052080           7.00307486       99996.12998843 
                                                    3.99982965           2.99988385           8.00050879           0.00000103 
                                                 3999.84513206        6632.74319101      995483.45650841           1.13082955         
                                          


                                          Замечу, что вариант по теореме здорово «нашумел» погрешностями.
                                            0
                                            На мой взгляд, накопление ошибки здесь берётся из делений и умножений при вычислении коэффициентов характеристического уравнения. В обычном способе деление идёт только в финале.
                                              +1
                                              Тут может существенную погрешность вносить процесс возведения матрицы в степень. Если для этой матрицы не соблюден критерий устойчивости Ляпунова (действительная часть хотя бы одного собственного числа положительна), малые возмущения, которые сопровождают каждую арифметическую операцию, «раскачают» нашу линейную систему так, что в итоге мы получим матрицу, далекую от реальной матрицы A^n.

                                              В методе же с союзной матрицей имеет место быть определенная симметрия. Очень возможно, что погрешность, выскочившая при получении одного минора в процессе вычисления вычитается из точно такой же погрешности, выскочившей при вычислении другого минора. Происходит своеобразная «стабилизация», компенсация погрешностей. Но это нужно скрупулезно доказывать.
                                          0
                                          Я небольшое отступление сделаю. Из применения теоремы Гамильтона-Кэли к обращению матриц получается интересная картина. Когда мы говорим о существовании обратной матрицы, предполагаем, что определитель матрицы A отличен от нуля (или все собственные числа этой матрицы должны быть ненулевыми). А тут получается, что для существования обратной матрицы нужно, чтобы коэффициент характеристического уравнения матрицы A^n был отличен от нуля.
                                            0
                                            Матрицы A, я описался.
                                              +2
                                              Сильно подозреваю, что можно доказать: этот коэффициент на самом деле и есть определитель, возможно умноженный на константу.
                                                +1
                                                Сейчас я прикинул, это и есть определитель, взятый с противоположным знаком (при записи и вычислениях я использовал перед ним знак "-"), поскольку является свободным членом в этом уравнении, а свободный член по теореме Виета равен произведению корней полинома, произведение характеристических чисел матрицы равно определителю. Ещё один способ его вычисления.
                          0
                          В C++14 в constexpr можно делать конструкторы и нормальные функции (в 11 только однострочные).
                          И это очень круто, если вместо const можно написать constexpr, тогда можно вовсе не парися о оптимизациях преде компиляцией, считая самому матрицы, и быть увереным что компилятор сам сложит выражения типа единичная матрица*перенос*поворот*инверс(перенос) в одну матрицу еще на компиляции.
                            +1
                            Нарисованное перед функцией constexpr вовсе не гарантирует, что лентяй-компилятор вычислит ее значение во время компиляции. Гарантию дает только использование std::integral_constant:
                            cout<<std::integral_constant<int, g(0)>::value;
                            
                              0
                              При чем тут вообще это? constexpr показывает что счисление возможно провести на этапе компиляции. Но это так же значит что функции помеченые constexpr можно вызывать и в коде не с constexpr выражения.
                              Библиотека линейной алгебры которая обязана все считать на компиляции мягко говоря, бесполезна ведь.
                              Да и вообще, меньше слов, больше кода — ideone.com/wRz8nr
                              Этот мини класс можно использовать и для генерации компайл-тайм констант и для любых рантайм данных.
                                –1
                                Очень интересно, особенно строка 15:
                                return std::move ( type(lhs.x + rhs.x, lhs.y + rhs.y) );
                                


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

                                Коротко: return std::move(x) — бессмысленная избыточность. Достаточно просто return x;

                                Это прямо описано в стандарте (12.8/32).
                                Обсуждение на stackoverflow
                                  0
                                  Снова не по теме. Это кусок кода который я написал скопировал из другого места, где std::move нужен был. Он демострирует НЕ std::move, который здесь вообще ничего не делает, и я это знаю.
                                  Смысл первого моего коментария в том что было бы интересно сделать полноценную компайл-тайм адекватную библиотеку матриц для задачи туториала.
                                  0
                                  constexpr показывает что счисление возможно провести на этапе компиляции.

                                  Так вам вроде об этом и сказали — только возможно, а не обязательно. Впрочем, я стандарт не читал, вдруг там гарантируется, что constexpr-функции от constexpr-выражений должны быть посчитаны на этапе компиляции? Это было бы здорово! Был бы рад, если кто-нибудь просветит по этому вопросу.
                                    0
                                    Не гарантирует этого стандарт, можно написать банальный тест и проверить. Только если результат пишется в constexpr — переменную, магия сработает. Собственно это и делает std::integral_constant()
                                      0
                                      Ну, тест это не гарантия, от багов никто не застрахован. Да и разработчики компиляторов любят вносить свои улучшения/отклонения от стандарта.

                                      Только если результат пишется в constexpr — переменную, магия сработает.

                                      Почему вот, кстати? Вроде должно быть очевидно, что это избыточное условие.
                                        –1
                                        Гипотетически(!), можно предположить, что компилятор в таком случае резонно поинтересуется:
                                        -А результат вычисления constexpr-функции я куда складировать буду? Вот переменную мне объявите, я в нее все сложу и в секцию статических данных выполняемого файла сохраню.

                                        То есть здесь на лицо использование уже имеющегося в языке функционала для размещения статических данных.
                                          0
                                          если просто, то constexpr бывает в двух ипостасях:
                                          — на переменных
                                          — на функциях
                                          Если в выражении используются только constexpr переменные, с constexpr фуркциями и литералы, то это выражение может быть вычислено в компайл тайме.
                                          Любое такое выражение можно исползьвать там где нужна константа компиляции.
                                          Так же функции constexpr можно использовать с динамическими данными, т.е. это никак не деградирует рантайм полезность библиотеки.
                                          0
                                          std::integral_constant это вообще тип, если что… Я окончательно перестал понимать что вы хотите сказать. Вам хочется что бы обязательно счисления проводились в компайл тайме? Ну тогда да, вы можете это зафорсить обвернув в этот тип, но все использованные функции Обязаны быть constexpr. Это не или/или.
                                          0
                                          Уф… Блин, походу тут какое-то недоразумение. constexpr показывает что функцию/объект можно использовать в качестве компайл тайм констант. Учитывая специфики современных компляторов — они сделают все что возможно в компайл тайме. То есть когда вы пишите int x = 5+4; вы ведь правда считаете что компилятор постчитает выражение 5+4 и вам не нужно прибавлять их самому? Когда делаете объекты и функции constexpr, это показывает компилятору так же что ваш объект можно так оптимизировать. Конечно, если полностью отключить оптимизацию он этого не сделает.

                                          Да, еще плюс в том что функции помеченые constexpr так же работают и на динамических данных. Т.е. в рантайме те же функции будут так же работать.
                                            0
                                            GCC не лентяйничает с constexpr только в режиме жесткой оптимизации (-O3).
                                            На -O2 он функцию только заинлайнил:
                                            #include <iostream>
                                            using namespace std;
                                            constexpr unsigned int fact(unsigned int in)
                                            {
                                              return  in==0
                                                      ?
                                                         1
                                                      :
                                                         in*fact(in-1);
                                            }
                                            int main() 
                                            {
                                                cout << fact(8) <<'\n';
                                                return 0;
                                            }
                                            

                                            Компилируем:
                                            g++ test.cpp -std=c++11 -O2 -S -o out.s
                                            

                                            Смотрим листинг:
                                            main:
                                            .LFB3682:
                                            	.cfi_startproc
                                            	movl	$1, %esi
                                            	movl	$8, %eax
                                            	.p2align 4,,10
                                            	.p2align 3
                                            .L3:
                                            	imull	%eax, %esi		;заинлайненное вычисление факториала.
                                            	subl	$1, %eax		;GCC - умница, развернул рекурсию
                                            	jne	.L3													;
                                            	subq	$24, %rsp
                                            	.cfi_def_cfa_offset 32
                                            	movl	$_ZSt4cout, %edi
                                            	call	_ZNSo9_M_insertImEERSoT_
                                            	leaq	15(%rsp), %rsi
                                            

                                            Компилируем с более жесткой оптимизацией:
                                            g++ test.cpp -std=c++11 -O3 -S -o out.s
                                            

                                            .cfi_startproc
                                            	subq	$24, %rsp
                                            	.cfi_def_cfa_offset 32
                                            	movl	$40320, %esi              ;все честно, 40320 - это факториал 8
                                            	movl	$_ZSt4cout, %edi
                                            	call	_ZNSo9_M_insertImEERSoT_
                                            	leaq	15(%rsp), %rsi
                                            

                                            Посчитали факториал в maxima:
                                            Maxima 5.35.1 http://maxima.sourceforge.net
                                            using Lisp CLISP 2.49 (2010-07-07)
                                            Distributed under the GNU Public License. See the file COPYING.
                                            Dedicated to the memory of William Schelter.
                                            The function bug_report() provides bug reporting information.
                                            (%i1) 8!;
                                            (%o1)                                40320
                                            (%i2) 
                                            


                                            Напоминаю, что по умолчанию в большинстве IDE стоит -O2 для релиза
                                              0
                                              А что будет если написать
                                              constexpr auto a = fact(8);
                                              cout << a <<'\n';
                                              вместо
                                              cout << fact(8); <<'\n';
                                              ?
                                                0
                                                В любом режиме оптимизации считает в compile-time, как и положено по стандарту.

                                                Итого: Чтобы применить вашу идею, и чтобы это было переносимо между разными компиляторами и их режимами, придется все такие места внимательно оборачивать в std::integral_constant.

                                                Иначе получится, например, что при переключении IDE в Debug (-O0) будет жуткая посадка производительности.
                                                  –1
                                                  Общим, похоже у вас единственная цель — доказать свою правоту, при чем не понятно в каком вопросе. Я не собираюсьразводить споров ради споров. У меня все работает при -O1. Дебаг — он для дебага. Он вообще не разворачивает ничего. Попробуйте в дебаге обвернуть что-то в std::integral_constant и сильно удивитесь.
                                                    +2
                                                    У меня цель — сделать корректно, надежно и переносимо.

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

                                                    Я замечаю, что constexpr вовсе не обязывает компилятор что-то там считать, ему нужно явно указать, что результат пишется в constexpr.

                                                    Кстати, в режиме отладки,
                                                    g++ test.cpp -std=c++11 -O0 -g -S -o out.s
                                                    

                                                    выражение обернутое в std::integral_constant
                                                    cout << std::integral_constant<unsigned int,fact(8)>::value;
                                                    

                                                    Прекрасно сосчиталось во время компиляции:
                                                    main:
                                                    .LFB3639:
                                                    	.file 1 "test.cpp"
                                                    	.loc 1 14 0
                                                    	.cfi_startproc
                                                    	pushq	%rbp
                                                    	.cfi_def_cfa_offset 16
                                                    	.cfi_offset 6, -16
                                                    	movq	%rsp, %rbp
                                                    	.cfi_def_cfa_register 6
                                                    	.loc 1 16 0
                                                    	movl	$40320, %esi                ;ага, вот эта константа
                                                    	movl	$_ZSt4cout, %edi
                                                    	call	_ZNSolsEj
                                                    	movl	$10, %esi
                                                    	movq	%rax, %rdi
                                                      0
                                                      1) А как выглядит дебаг для такого кода? Мы ведь правда хотим дебажить наши функции и только для этого делаем билды с -O0.

                                                      2) Кстате, заметте, что что бы использовать fact в integral_constant вам необходимо сделать его constexpr. Именно об этом я и говорил. Уже 5 раз. Что если сделать функцию constexpr — вы ее можете использовать как копайл константу, например в типе integral_constant. Но при этом это не ограничивает ее использвание. Мы так же может использовать fact от пользовательских данных. Пользовательские данные никак не обернуть в integral_constant, потому что это противоречит смыслу integral_constant.
                                                      Я вообще не понимаю что вы хотите сказать? Что обязательно все обворачивать в integral_constant? Спасибо, за меня это сделает даже О1 оптимизация.
                                                        +1
                                                        Вы правильно заметили — отлаживать код, который выполняется в compile-time невозможно. Также, как невозможно отлаживать, например, раскрутку рекурсии на шаблонах.

                                                        Я хочу сказать и показать, что для того, чтобы значение constexpr-функции гарантированно вычислилось в compile-time, его необходимо сохранять в статическую константу. Если этого не сделать, вычисление на этапе компиляции будет выполняться только в режиме агрессивной оптимизации -O3, а не -O1, как вы утверждаете:
                                                        Спасибо, за меня это сделает даже О1 оптимизация
                                                          0
                                                          Я это всегда знал ведь ;) Мы просто непоняли друг друга сначала.
                                                        0
                                                        Думаю мы поняли друг друга. Общим, все что я хочу это что бы вот что было:
                                                        constexpr int det = matrix([1,4,9],[2,9,7],[4,3,8]).det();
                                                        cout << std::integral_constant<int, det >::value;
                                                        ну или как-то так.
                                                          +2
                                                          После строки:
                                                          constexpr int det=...
                                                          

                                                          оборачивание det в std::integral_constant бессмысленно — запись в constexpr-переменную гарантирует вычисление в compile-time.

                                                          Я собираюсь сделать весь рендер во время компиляции, не только матричную алгебру.
                                                            0
                                                            Вобще смысл писать «constexpr int det» в том что можно просто удалить слово constexpr и тогда можно дебаггером зайти внутрь выражения. integral_constant действительно лишний семантически, он для выразительности что это точно просто константа в бинарнике.
                                                            Вообще очень клевая идея сделать рендерер в компиляторе. А как обойдетесь io?
                                                              0
                                                              Входные данные поместить в строку, выходные — либо в кортеж, либо в массив.
                                                      +2
                                                      В любом режиме оптимизации считает в compile-time, как и положено по стандарту.
                                                      Честно говоря, я не уверен, что стандарт где-то обязывает так делать. В конце концов, у нас есть правило «as-if», а поскольку программы с вычислением выражения в compile-time и в runtime имеют одно и то же наблюдаемое поведение, то они обе корректны. На практике, конечно, понятно, что когда мы используем выражение как параметр шаблона, то компилятор его посчитает во время компиляции, просто чтобы проверить, что не будет ошибки компиляции, дать имя этому шаблону и т.д. Тут было лучше какую-нибудь опцию компилятора иметь, которая говорила бы ему «считай все constexpr в compile-time», но для gcc в документации я такой не нашел.
                                                        0
                                                        Здесь можно за уши притянуть тот факт, что значение в constexpr должно быть доступно во время компиляции.

                                                        Но формально, если не заставлять компилятор силком (используя std::integral_constant, для которого это значение и будет параметром шаблона), компилятор может и слентяйничать.

                                                        Итого: с «положено по стандарту» я действительно погорячился.
                                                          +1
                                                          Мне кажется, с переменными та же самая история, что и с параметром шаблона. Во-первых, компилятор должен проверить, что значение, которым она инициализируется, вычислимо в compile-time, и если нет, то выдать ошибку компиляции. Причем для этой проверки недостаточно просто проверить, что функция и её аргумент constexpr, например:
                                                          int non_constexpr() { return 1; }
                                                          
                                                          constexpr int f(bool b)
                                                          {
                                                              return b ? 0 : non_constexpr();
                                                          }
                                                          
                                                          constexpr int x = f(true); // OK
                                                          constexpr int y = f(false); // Compilation error
                                                          

                                                          Во-вторых, переменная может быть использовано в контексте, требующем знания её значения в момент компиляции (параметр шаблона, размерность массива, static_assert()), поэтому, возможно, компилятору тупо проще вычислить значение сразу, чем проводить анализ, а будет ли оно нам нужно, или же сохранять какие-то данные для ленивого вычисления.
                                                    –1
                                                    Общим, у меня все работает на -O1 — goo.gl/XVhbOJ
                                                    Я ж не говорю что оно работает просто так. Constexpr тоже готовить нужно уметь. Вот для сравнения ваш код goo.gl/bR14sW

                                        Only users with full accounts can post comments. Log in, please.