Dingu'yddie

This commit is contained in:
william 2024-02-13 21:38:14 -05:00
parent 2a412ddb13
commit fcaa08a89a
6 changed files with 370 additions and 139 deletions

View File

@ -26,85 +26,114 @@
namespace gti320 {
// Deux types de vecteurs 3D considérés ici
typedef Vector<double, 3> Vector3d;
typedef Vector<float, 3> Vector3f;
// Deux types de vecteurs 3D considérés ici
typedef Vector<double, 3> Vector3d;
typedef Vector<float, 3> Vector3f;
// Dans le cadre de ce projet, nous considérons seulement deux
// cas :
//
// - les rotations
// - les translations
//
// Deux types de matrices en coordonnées homogèes :
typedef Matrix<double, 4, 4, ColumnStorage> Matrix4d;
typedef Matrix<float, 4, 4, ColumnStorage> Matrix4f;
//
// Deux types de matrices pour les rotations
typedef Matrix<double, 3, 3, ColumnStorage> Matrix3d;
typedef Matrix<float, 3, 3, ColumnStorage> Matrix3f;
// Dans le cadre de ce projet, nous considérons seulement deux
// cas :
//
// - les rotations
// - les translations
//
// Deux types de matrices en coordonnées homogèes :
typedef Matrix<double, 4, 4, ColumnStorage> Matrix4d;
typedef Matrix<float, 4, 4, ColumnStorage> Matrix4f;
//
// Deux types de matrices pour les rotations
typedef Matrix<double, 3, 3, ColumnStorage> Matrix3d;
typedef Matrix<float, 3, 3, ColumnStorage> Matrix3f;
/**
* Initialise et retourne la matrice identité
*/
template<>
inline void Matrix4d::setIdentity()
{
// TODO affecter la valeur 0.0 partout, sauf sur la diagonale principale où c'est 1.0.
// Note: ceci est une redéfinition d'une fonction membre!
}
/**
* Initialise et retourne la matrice identité
*/
template<>
inline void Matrix4d::setIdentity() {
this->setZero();
/**
* Calcul de la matrice inverse, SPÉCIALISÉ pour le cas d'une matrice de
* transformation en coordonnées homogènes.
*
* TODO (vous pouvez supposer qu'il s'agit d'une matrice de transformation
* en coordonnées homogènes)
*/
template<>
inline Matrix4d Matrix4d::inverse() const
{
// TODO : implémenter
return Matrix4d(); // Pas bon, à changer
}
for (int i = 0; i < 4; i++) {
(*this)(i, i) = 1;
}
}
/**
* Calcul de la matrice inverse, SPÉCIALISÉ pour le cas d'une matrice de rotation.
*
* (vous pouvez supposer qu'il s'agit d'une matrice de rotation)
*/
template<>
inline Matrix3d Matrix3d::inverse() const
{
// TODO : implémenter
return Matrix3d();
}
/**
* Calcul de la matrice inverse, SPÉCIALISÉ pour le cas d'une matrice de
* transformation en coordonnées homogènes.
*/
template<>
inline Matrix4d Matrix4d::inverse() const {
auto transposed_rotation = this->block(0, 0, 3, 3).transpose<double, 3, 3, ColumnStorage>();
auto translation = Vector<double, 3>();
for (int i = 0; i < 3; i++) {
translation(i) = (*this)(i, 3);
}
auto result = Matrix4d();
result(3, 3) = 1;
// Applique la matrice de rotation
result.block(0, 0, 3, 3) = transposed_rotation;
// Applique le vecteur de translation
auto inverted_translation = (-1.0 * transposed_rotation) * translation;
auto trans_block = result.block(0, 3, 3, 1);
trans_block(0, 0) = inverted_translation(0);
trans_block(1, 0) = inverted_translation(1);
trans_block(2, 0) = inverted_translation(2);
return result;
}
/**
* Calcul de la matrice inverse, SPÉCIALISÉ pour le cas d'une matrice de rotation.
*
* (vous pouvez supposer qu'il s'agit d'une matrice de rotation)
*/
template<>
inline Matrix3d Matrix3d::inverse() const {
return this->transpose<double, 3, 3, ColumnStorage>();
}
/**
* Multiplication d'une matrice 4x4 avec un vecteur 3D la matrice
* représente une transformation en coordonnées homogène.
*/
template<typename _Scalar>
Vector<_Scalar, 3> operator*(const Matrix<_Scalar, 4, 4, ColumnStorage> &A, const Vector<_Scalar, 3> &v) {
auto homogeneous = Vector<_Scalar, 4>(v);
homogeneous(3) = 1;
Vector<_Scalar, 4> transformed = A * homogeneous;
return Vector<_Scalar, 3>(transformed);
}
/**
* Multiplication d'une matrice 4x4 avec un vecteur 3D la matrice
* représente une transformation en coordonnées homogène.
*/
template <typename _Scalar>
Vector<_Scalar, 3> operator*(const Matrix<_Scalar, 4, 4, ColumnStorage>& A, const Vector<_Scalar, 3>& v)
{
// TODO : implémenter
return Vector<_Scalar, 3>(); // pas bon, à changer
}
/**
* Initialise et retourne la matrice de rotation définie par les angles
* d'Euler XYZ exprimés en radians.
*
* La matrice doit correspondre au produit : Rz*Ry*Rx.
*/
template<typename _Scalar>
static Matrix<_Scalar, 3, 3> makeRotation(_Scalar x, _Scalar y, _Scalar z) {
double sin_theta = std::sin(x);
double cos_theta = std::cos(x);
double sin_phi = std::sin(y);
double cos_phi = std::cos(y);
double sin_psi = std::sin(z);
double cos_psi = std::cos(z);
auto m = Matrix<_Scalar, 3, 3>();
m(0, 0) = cos_phi * cos_psi;
m(1, 0) = cos_phi * sin_psi;
m(2, 0) = -sin_phi;
m(0, 1) = sin_theta * sin_phi * cos_psi - cos_theta * sin_psi;
m(1, 1) = sin_theta * sin_phi * sin_psi + cos_theta * cos_psi;
m(2, 1) = sin_theta * cos_phi;
m(0, 2) = cos_theta * sin_phi * cos_psi + sin_theta * sin_psi;
m(1, 2) = cos_theta * sin_phi * sin_psi - sin_theta * cos_psi;
m(2, 2) = cos_theta * cos_phi;
/**
* Initialise et retourne la matrice de rotation définie par les angles
* d'Euler XYZ exprimés en radians.
*
* La matrice doit correspondre au produit : Rz*Ry*Rx.
*/
template<typename _Scalar>
static Matrix<_Scalar, 3, 3> makeRotation(_Scalar x, _Scalar y, _Scalar z)
{
// TODO : implémenter
return Matrix<_Scalar, 3, 3>(); // pas bon, à changer
}
return m;
}
}

View File

@ -12,6 +12,7 @@
*/
#include <complex>
#include <iostream>
#include "MatrixBase.h"
namespace gti320 {

View File

@ -37,9 +37,12 @@ namespace gti320
/**
* Constructeur de copie
*/
MatrixBase(const MatrixBase& other) : m_storage(other.m_storage) { }
MatrixBase(const MatrixBase& other) : m_storage(other.m_storage) {
}
explicit MatrixBase(int _rows, int _cols) : m_storage() { }
explicit MatrixBase(int _rows, int _cols) : m_storage() {
setZero();
}
/**
* Destructeur
@ -119,7 +122,9 @@ namespace gti320
*/
MatrixBase() : m_storage(), m_rows(0) { }
explicit MatrixBase(int _rows, int _cols) : m_storage(_rows* _Cols), m_rows(_rows) { }
explicit MatrixBase(int _rows, int _cols) : m_storage(_rows* _Cols), m_rows(_rows) {
setZero();
}
/**
* Constructeur de copie
@ -205,7 +210,9 @@ namespace gti320
*/
MatrixBase() : m_storage() { }
explicit MatrixBase(int _rows, int _cols) : m_storage(_rows* _cols), m_cols(_cols) { }
explicit MatrixBase(int _rows, int _cols) : m_storage(_rows* _cols), m_cols(_cols) {
setZero();
}
/**
* Constructeur de copie
@ -291,7 +298,9 @@ namespace gti320
*/
MatrixBase() : m_storage(), m_rows(0), m_cols(0) { }
explicit MatrixBase(int _rows, int _cols) : m_storage(_rows* _cols), m_rows(_rows), m_cols(_cols) { }
explicit MatrixBase(int _rows, int _cols) : m_storage(_rows* _cols), m_rows(_rows), m_cols(_cols) {
setZero();
}
/**
* Constructeur de copie

View File

@ -127,7 +127,7 @@ namespace gti320 {
for (int col = 0; col < A.cols(); col++) {
for (int row = 0; row < A.rows(); row++) {
result(col, row) = A(col, row) + B(col, row);
result(row, col) = A(row, col) + B(row, col);
}
}
@ -206,8 +206,21 @@ namespace gti320 {
template<typename _Scalar, int _Rows, int _Cols>
Vector<_Scalar, _Rows>
operator*(const Matrix<_Scalar, _Rows, _Cols, RowStorage> &A, const Vector<_Scalar, _Cols> &v) {
// TODO : implémenter
return Vector<_Scalar, _Rows>();
assert(A.cols() == v.size());
auto result = Vector<_Scalar, _Rows>(A.rows());
for (int row = 0; row < A.rows(); row++) {
_Scalar dotP = 0;
for (int col = 0; col < A.cols(); col++) {
dotP += A(col, row) * v(col);
}
result(row) = dotP;
}
return result;
}
/**
@ -219,8 +232,19 @@ namespace gti320 {
template<typename _Scalar, int _Rows, int _Cols>
Vector<_Scalar, _Rows>
operator*(const Matrix<_Scalar, _Rows, _Cols, ColumnStorage> &A, const Vector<_Scalar, _Cols> &v) {
// TODO : implémenter
return Vector<_Scalar, _Rows>();
assert(A.cols() == v.size());
auto result = Vector<_Scalar, _Rows>(A.rows());
for (int col = 0; col < A.cols(); col++) {
_Scalar v_col = v(col);
for (int row = 0; row < A.rows(); row++) {
result(row) += A(row, col) * v_col;
}
}
return result;
}
/**
@ -228,8 +252,13 @@ namespace gti320 {
*/
template<typename _Scalar, int _Rows>
Vector<_Scalar, _Rows> operator*(const _Scalar &a, const Vector<_Scalar, _Rows> &v) {
// TODO : implémenter
return Vector<_Scalar, _Rows>();
auto result = Vector<_Scalar, _Rows>(v.rows());
for (int row = 0; row < v.rows(); row++) {
result(row) = a * v(row);
}
return result;
}
@ -238,8 +267,15 @@ namespace gti320 {
*/
template<typename _Scalar, int _RowsA, int _RowsB>
Vector<_Scalar, _RowsA> operator+(const Vector<_Scalar, _RowsA> &a, const Vector<_Scalar, _RowsB> &b) {
// TODO : implémenter
return Vector<_Scalar, _RowsA>();
assert(a.rows() == b.rows());
auto result = Vector<_Scalar, _RowsA>(a.rows());
for (int row = 0; row < a.rows(); row++) {
result(row) = a(row) + b(row);
}
return result;
}
/**
@ -247,7 +283,82 @@ namespace gti320 {
*/
template<typename _Scalar, int _RowsA, int _RowsB>
Vector<_Scalar, _RowsA> operator-(const Vector<_Scalar, _RowsA> &a, const Vector<_Scalar, _RowsB> &b) {
// TODO : implémenter
return Vector<_Scalar, _RowsA>();
assert(a.rows() == b.rows());
auto result = Vector<_Scalar, _RowsA>(a.rows());
for (int row = 0; row < a.rows(); row++) {
result(row) = a(row) - b(row);
}
return result;
}
template<typename _Scalar, int _Rows, int _Cols>
void print_matrix(const Matrix<_Scalar, _Rows, _Cols> &m) {
std::cout << "+----" << '\n';
std::cout.precision(2);
for (int row = 0; row < m.rows(); row++) {
for (int col = 0; col < m.cols(); col++) {
std::cout << m(row, col) << "," << "\t";
}
std::cout << "\n";
}
}
template<typename _Scalar>
_Scalar det(const Matrix<_Scalar, 1, 1> &m) {
return m(0, 0);
}
template<typename _Scalar>
_Scalar det(const Matrix<_Scalar, 2, 2> &m) {
_Scalar a = m(0, 0);
_Scalar b = m(0, 1);
_Scalar c = m(1, 0);
_Scalar d = m(1, 1);
return (a * d) - (b * c);
}
template<typename _Scalar, int _Rows, int _Cols>
double det(const Matrix<_Scalar, _Rows, _Cols> &m) {
assert(m.rows() == m.cols());
auto l = Matrix<double, _Rows, _Cols>();
auto u = Matrix<double, _Rows, _Cols>();
for (int j = 0; j < m.cols(); j++) {
for (int i = 0; i < m.rows(); i++) {
u(i, j) = m(i, j);
l(i, j) = m(i, j);
for (int k = 0; k < i; k++) {
u(i, j) -= l(i, k) * u(k, j);
}
for (int k = 0; k < j; k++) {
l(i, j) -= l(i, k) * u(k, j);
}
l(i, j) = 1 / u(j, j);
}
}
print_matrix(l);
print_matrix(u);
return det_tri(l) * det_tri(u);
}
template<typename _Scalar, int _Rows, int _Cols>
_Scalar det_tri(const Matrix<_Scalar, _Rows, _Cols> &m) {
int n = std::min(m.rows(), m.cols());
int det = 1;
for (int i = 0; i < n; i++) {
det *= m(i, i);
}
return det;
}
}

View File

@ -41,6 +41,23 @@ namespace gti320 {
*/
Vector(const Vector &other) : MatrixBase<_Scalar, _Rows, 1>(other) {}
/**
* Constructeur de copie avec une grandeur différente
* @tparam _OtherRows
* @param other
*/
template<int _OtherRows>
Vector(const Vector<_Scalar, _OtherRows> &other) : MatrixBase<_Scalar, _Rows, 1>() {
int n = std::min(other.size(), _Rows);
const _Scalar *src = other.data();
_Scalar *dest = this->m_storage.data();
memcpy(dest, src, n * sizeof(_Scalar));
if (n < _Rows) {
memset(&dest[n], 0, (_Rows - n) * sizeof(_Scalar));
}
}
/**
* Destructeur
*/

View File

@ -171,64 +171,84 @@ TEST(TestLabo1, Matrix4x4SizeTest) {
*/
TEST(TestLabo1, MatrixMatrixOperators) {
// Opérations arithmétiques avec matrices à taille dynamique
{
// Test : matrice identité
Matrix<double> A(6, 6);
A.setIdentity();
EXPECT_DOUBLE_EQ(A(0, 0), 1.0);
EXPECT_DOUBLE_EQ(A(1, 1), 1.0);
EXPECT_DOUBLE_EQ(A(2, 2), 1.0);
EXPECT_DOUBLE_EQ(A(3, 3), 1.0);
EXPECT_DOUBLE_EQ(A(4, 4), 1.0);
EXPECT_DOUBLE_EQ(A(5, 5), 1.0);
EXPECT_DOUBLE_EQ(A(0, 1), 0.0);
EXPECT_DOUBLE_EQ(A(1, 0), 0.0);
{
// Test : matrice identité
Matrix<double> A(6, 6);
A.setIdentity();
EXPECT_DOUBLE_EQ(A(0, 0), 1.0);
EXPECT_DOUBLE_EQ(A(1, 1), 1.0);
EXPECT_DOUBLE_EQ(A(2, 2), 1.0);
EXPECT_DOUBLE_EQ(A(3, 3), 1.0);
EXPECT_DOUBLE_EQ(A(4, 4), 1.0);
EXPECT_DOUBLE_EQ(A(5, 5), 1.0);
EXPECT_DOUBLE_EQ(A(0, 1), 0.0);
EXPECT_DOUBLE_EQ(A(1, 0), 0.0);
// Test : produit scalaire * matrice
const double alpha = 2.5;
Matrix<double> B = alpha * A;
EXPECT_DOUBLE_EQ(B(0, 0), alpha);
EXPECT_DOUBLE_EQ(B(1, 1), alpha);
EXPECT_DOUBLE_EQ(B(2, 2), alpha);
EXPECT_DOUBLE_EQ(B(3, 3), alpha);
EXPECT_DOUBLE_EQ(B(4, 4), alpha);
EXPECT_DOUBLE_EQ(B(5, 5), alpha);
EXPECT_DOUBLE_EQ(B(0, 1), 0.0);
EXPECT_DOUBLE_EQ(B(1, 0), 0.0);
// Test : produit scalaire * matrice
const double alpha = 2.5;
Matrix<double> B = alpha * A;
EXPECT_DOUBLE_EQ(B(0, 0), alpha);
EXPECT_DOUBLE_EQ(B(1, 1), alpha);
EXPECT_DOUBLE_EQ(B(2, 2), alpha);
EXPECT_DOUBLE_EQ(B(3, 3), alpha);
EXPECT_DOUBLE_EQ(B(4, 4), alpha);
EXPECT_DOUBLE_EQ(B(5, 5), alpha);
EXPECT_DOUBLE_EQ(B(0, 1), 0.0);
EXPECT_DOUBLE_EQ(B(1, 0), 0.0);
// Test : produit matrice * matrice
Matrix<double> C = A * B;
EXPECT_DOUBLE_EQ(C(0, 0), A(0, 0) * B(0, 0));
EXPECT_DOUBLE_EQ(C(1, 1), A(1, 1) * B(1, 1));
EXPECT_DOUBLE_EQ(C(2, 2), A(2, 2) * B(2, 2));
EXPECT_DOUBLE_EQ(C(3, 3), A(3, 3) * B(3, 3));
EXPECT_DOUBLE_EQ(C(4, 4), A(4, 4) * B(4, 4));
EXPECT_DOUBLE_EQ(C(5, 5), A(5, 5) * B(5, 5));
EXPECT_DOUBLE_EQ(C(0, 1), 0.0);
EXPECT_DOUBLE_EQ(C(2, 3), 0.0);
// Test : produit matrice * matrice
Matrix<double> C = A * B;
EXPECT_DOUBLE_EQ(C(0, 0), A(0, 0) * B(0, 0));
EXPECT_DOUBLE_EQ(C(1, 1), A(1, 1) * B(1, 1));
EXPECT_DOUBLE_EQ(C(2, 2), A(2, 2) * B(2, 2));
EXPECT_DOUBLE_EQ(C(3, 3), A(3, 3) * B(3, 3));
EXPECT_DOUBLE_EQ(C(4, 4), A(4, 4) * B(4, 4));
EXPECT_DOUBLE_EQ(C(5, 5), A(5, 5) * B(5, 5));
EXPECT_DOUBLE_EQ(C(0, 1), 0.0);
EXPECT_DOUBLE_EQ(C(2, 3), 0.0);
// Test : addition matrice * matrice
Matrix<double> A_plus_B = A + B;
EXPECT_DOUBLE_EQ(A_plus_B(0, 0), A(0, 0) + B(0, 0));
EXPECT_DOUBLE_EQ(A_plus_B(1, 1), A(1, 1) + B(1, 1));
EXPECT_DOUBLE_EQ(A_plus_B(2, 2), A(2, 2) + B(2, 2));
EXPECT_DOUBLE_EQ(A_plus_B(3, 3), A(3, 3) + B(3, 3));
EXPECT_DOUBLE_EQ(A_plus_B(4, 4), A(4, 4) + B(4, 4));
EXPECT_DOUBLE_EQ(A_plus_B(5, 5), A(5, 5) + B(5, 5));
EXPECT_DOUBLE_EQ(A_plus_B(0, 1), 0.0);
EXPECT_DOUBLE_EQ(A_plus_B(2, 3), 0.0);
}
// Test : addition matrice + matrice
Matrix<double> A_plus_B = A + B;
EXPECT_DOUBLE_EQ(A_plus_B(0, 0), A(0, 0) + B(0, 0));
EXPECT_DOUBLE_EQ(A_plus_B(1, 1), A(1, 1) + B(1, 1));
EXPECT_DOUBLE_EQ(A_plus_B(2, 2), A(2, 2) + B(2, 2));
EXPECT_DOUBLE_EQ(A_plus_B(3, 3), A(3, 3) + B(3, 3));
EXPECT_DOUBLE_EQ(A_plus_B(4, 4), A(4, 4) + B(4, 4));
EXPECT_DOUBLE_EQ(A_plus_B(5, 5), A(5, 5) + B(5, 5));
EXPECT_DOUBLE_EQ(A_plus_B(0, 1), 0.0);
EXPECT_DOUBLE_EQ(A_plus_B(2, 3), 0.0);
}
// Opérations arithmétique avec matrices à stockage par lignes et par
// colonnes.
{
// Création d'un matrice à stockage par lignes
Matrix<double, Dynamic, Dynamic, RowStorage> A(5, 5);
A(0, 0) = 0.8147; A(0, 1) = 0.0975; A(0, 2) = 0.1576; A(0, 3) = 0.1419; A(0, 4) = 0.6557;
A(1, 0) = 0.9058; A(1, 1) = 0.2785; A(1, 2) = 0.9706; A(1, 3) = 0.4218; A(1, 4) = 0.0357;
A(2, 0) = 0.1270; A(2, 1) = 0.5469; A(2, 2) = 0.9572; A(2, 3) = 0.9157; A(2, 4) = 0.8491;
A(3, 0) = 0.9134; A(3, 1) = 0.9575; A(3, 2) = 0.4854; A(3, 3) = 0.7922; A(3, 4) = 0.9340;
A(4, 0) = 0.6324; A(4, 1) = 0.9649; A(4, 2) = 0.8003; A(4, 3) = 0.9595; A(4, 4) = 0.6787;
A(0, 0) = 0.8147;
A(0, 1) = 0.0975;
A(0, 2) = 0.1576;
A(0, 3) = 0.1419;
A(0, 4) = 0.6557;
A(1, 0) = 0.9058;
A(1, 1) = 0.2785;
A(1, 2) = 0.9706;
A(1, 3) = 0.4218;
A(1, 4) = 0.0357;
A(2, 0) = 0.1270;
A(2, 1) = 0.5469;
A(2, 2) = 0.9572;
A(2, 3) = 0.9157;
A(2, 4) = 0.8491;
A(3, 0) = 0.9134;
A(3, 1) = 0.9575;
A(3, 2) = 0.4854;
A(3, 3) = 0.7922;
A(3, 4) = 0.9340;
A(4, 0) = 0.6324;
A(4, 1) = 0.9649;
A(4, 2) = 0.8003;
A(4, 3) = 0.9595;
A(4, 4) = 0.6787;
// Test : transposée (le résultat est une matrice à stockage par
// colonnes)
@ -486,8 +506,6 @@ TEST(TestLabo1, Math3D) {
EXPECT_DOUBLE_EQ(Rinv(0, 0), RT(0, 0));
EXPECT_DOUBLE_EQ(Rinv(1, 1), RT(1, 1));
EXPECT_DOUBLE_EQ(Rinv(0, 2), RT(0, 2));
}
/**
@ -630,6 +648,52 @@ TEST(TestLabo1, Supplementaires) {
V1 = V2;
EXPECT_EQ(V1.size(), V2.size());
EXPECT_EQ(V1(12), V2(12));
// === Copie de vecteurs de grandeur différente ===
// Test 11: Petit vers grand
Vector<int, 5> V3;
for (int i = 0; i < V3.size(); i++) {
V3(i) = i;
}
Vector<int, 7> V4(V3);
for (int i = 0; i < V3.size(); i++) {
EXPECT_EQ(V4(i), V3(i));
}
for (int i = V3.size(); i < V4.size(); i++) {
EXPECT_EQ(V4(i), 0);
}
// Test 12: Grand vers petit
Vector<int, 5> V5(V4);
for (int i = 0; i < V5.size(); i++) {
EXPECT_EQ(V5(i), V4(i));
}
// === Déterminants ===
Matrix<int, 1, 1> MD1x1;
MD1x1(0, 0) = 1;
Matrix<int, 2, 2> MD2x2;
MD2x2(0, 0) = 1; MD2x2(0, 1) = 2;
MD2x2(1, 0) = 3; MD2x2(1, 1) = 4;
Matrix<int, 3, 3> MD3x3;
MD3x3(0, 0) = 1; MD3x3(0, 1) = 2; MD3x3(0, 2) = 3;
MD3x3(1, 0) = 4; MD3x3(1, 1) = 5; MD3x3(1, 2) = 6;
MD3x3(2, 0) = 7; MD3x3(2, 1) = 9; MD3x3(2, 2) = 0;
// Test 13: 1x1
int det1 = det(MD1x1);
EXPECT_EQ(det1, MD1x1(0, 0));
// Test 14: 2x2
int det2 = det(MD2x2);
EXPECT_EQ(det2, -2);
// Test 15: 3x3
int det3 = det(MD3x3);
EXPECT_EQ(det3, 33);
}
int main(int argc, char **argv) {