nerdcruft/libcruft-util
1
0
Fork 0
Code Issues Releases Activity

maths: remove deprecated maths objects

Browse Source
This commit is contained in:
Danny Robson 2015-06-05 16:07:03 +10:00
parent 62f97f0ec6
commit e8d3cf8eb1
6 changed files with 0 additions and 927 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
Makefile.am
View File

@ -121,10 +121,6 @@ UTIL_FILES = \
maths.cpp \
maths.hpp \
maths.ipp \
maths/matrix.cpp \
maths/matrix.hpp \
maths/vector.cpp \
maths/vector.hpp \
matrix.cpp \
matrix.hpp \
matrix.ipp \
@ -322,7 +318,6 @@ TEST_BIN = \
test/json_types \
test/ray \
test/maths \
test/maths_matrix \
test/matrix \
test/md2 \
test/md4 \

530
maths/matrix.cpp
View File

@ -1,530 +0,0 @@
/*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
* Copyright 2010 Danny Robson <danny@nerdcruft.net>
*/
#include "matrix.hpp"
#include "debug.hpp"
#include "range.hpp"
#include "maths.hpp"
#include <algorithm>
using namespace util;
using namespace maths;
matrix::matrix (size_t _rows, size_t _columns):
m_rows (_rows),
m_columns (_columns),
m_data (new double[_rows * _columns])
{
}
matrix::matrix (size_t _rows,
size_t _columns,
const std::initializer_list <double> &_data):
m_rows (_rows),
m_columns (_columns)
{
if (size () != _data.size ())
throw std::runtime_error ("element and initializer size differs");
CHECK_EQ (m_rows * m_columns, _data.size());
m_data.reset (new double[size ()]);
std::copy (_data.begin (), _data.end (), m_data.get ());
}
matrix::matrix (const std::initializer_list <vector> &rhs):
m_rows (rhs.size ()),
m_columns (rhs.begin()->size ()),
m_data (new double[m_rows * m_columns])
{
double *row_cursor = m_data.get ();
for (auto i = rhs.begin (); i != rhs.end (); ++i) {
CHECK (i->size () == m_columns);
std::copy (i->data (), i->data () + i->size (), row_cursor);
row_cursor += m_columns;
}
}
matrix::matrix (const matrix &rhs):
m_rows (rhs.m_rows),
m_columns (rhs.m_columns)
{
m_data.reset (new double [m_rows * m_columns]);
std::copy (rhs.m_data.get (), rhs.m_data.get () + m_rows * m_columns, m_data.get ());
}
matrix::matrix (matrix &&rhs):
m_rows (rhs.m_rows),
m_columns (rhs.m_columns),
m_data (std::move (rhs.m_data))
{
}
matrix::~matrix()
{ ; }
void
matrix::sanity (void) const {
CHECK (m_rows > 0);
CHECK (m_columns > 0);
CHECK (m_data != nullptr);
}
const double *
matrix::operator [] (unsigned int row) const {
CHECK_LT (row, m_rows);
return m_data.get () + row * m_columns;
}
double *
matrix::operator [] (unsigned int row) {
CHECK_LT (row, m_rows);
return m_data.get () + row * m_columns;
}
const double *
matrix::data (void) const
{ return m_data.get (); }
matrix&
matrix::operator =(const matrix& rhs) {
if (size () != rhs.size ()) {
m_data.reset (new double [rhs.rows () * rhs.columns ()]);
}
m_rows = rhs.m_rows;
m_columns = rhs.m_columns;
std::copy (rhs.m_data.get (), rhs.m_data.get () + m_rows * m_columns, m_data.get ());
return *this;
}
matrix
matrix::operator * (double scalar) const {
matrix val (*this);
for (unsigned int i = 0; i < m_rows; ++i)
for (unsigned int j = 0; j < m_columns; ++j)
val[i][j] *= scalar;
return val;
}
matrix&
matrix::operator *=(double scalar) {
for (unsigned int i = 0; i < m_rows; ++i)
for (unsigned int j = 0; j < m_columns; ++j)
(*this)[i][j] *= scalar;
return *this;
}
matrix&
matrix::operator /= (double scalar)
{ return (*this) *= (1.0 / scalar); }
matrix
matrix::operator + (double scalar) const {
matrix val (*this);
for (unsigned int i = 0; i < m_rows; ++i)
for (unsigned int j = 0; j < m_columns; ++j)
val[i][j] += scalar;
return val;
}
matrix&
matrix::operator +=(double scalar) {
for (unsigned int i = 0; i < m_rows; ++i)
for (unsigned int j = 0; j < m_columns; ++j)
(*this)[i][j] += scalar;
return *this;
}
matrix
matrix::operator * (const matrix& rhs) const {
if (m_columns != rhs.rows ())
throw std::invalid_argument ("matrices size mismatch in multiplication");
matrix val (matrix::zeroes (m_rows, rhs.columns ()));
for (unsigned int i = 0; i < m_rows; ++i)
for (unsigned int j = 0; j < rhs.columns (); ++j)
for (unsigned int k = 0; k < m_columns; ++k)
val[i][j] += (*this)[i][k] * rhs[k][j];
return val;
}
matrix&
matrix::operator *=(const matrix& rhs)
{ return *this = *this * rhs; }
bool
matrix::operator ==(const matrix& rhs) const {
if (rhs.rows () != rows () ||
rhs.columns () != columns ())
return false;
return std::equal (m_data.get (), m_data.get () + size (), rhs.data ());
}
//matrix transpose (void) const { ; }
size_t
matrix::rows (void) const
{ return m_rows; }
size_t
matrix::columns (void) const
{ return m_columns; }
size_t
matrix::size (void) const
{ return rows () * columns (); }
bool
matrix::is_square (void) const
{ return m_rows == m_columns; }
bool
matrix::is_magic (void) const {
if (!is_square ())
return false;
unsigned int expected = m_rows * (m_rows * m_rows + 1) / 2;
range<double> numbers (1, m_rows * m_rows);
for (unsigned int i = 0; i < m_rows; ++i) {
unsigned int sum1 = 0, sum2 = 0;
for (unsigned int j = 0; j < m_columns; ++j) {
double a = (*this)[i][j],
b = (*this)[j][i];
if (!exactly_equal (static_cast<uintmax_t> (a), a) ||
!exactly_equal (static_cast<uintmax_t> (b), b))
return false;
if (!numbers.contains (a) ||
!numbers.contains (b))
return false;
sum1 += static_cast<unsigned> (a);
sum2 += static_cast<unsigned> (b);
}
if (sum1 != expected || sum2 != expected)
return false;
}
return true;
}
bool
matrix::is_homogeneous (void) const {
if (m_rows != m_columns)
return false;
// Check the final row is all zeroes
for (unsigned int i = 0; i < m_columns - 1; ++i) {
if (!almost_equal ((*this)[m_rows - 1][i], 0.))
return false;
}
// Except for the last element, which has to be one
return almost_equal ((*this)[m_rows - 1][m_columns - 1], 1.);
}
double
matrix::determinant (void) const {
if (m_rows != m_columns)
not_implemented ();
switch (m_rows) {
case 2: return determinant2x2 ();
case 3: return determinant3x3 ();
case 4: return determinant4x4 ();
}
not_implemented ();
}
// With matrix A = [ a, b ]
// [ c, d ]
//
// det (A) = ad - bc
double
matrix::determinant2x2 (void) const {
CHECK_EQ (m_rows, 2);
CHECK_EQ (m_columns, 2);
return (*this)[0][0] * (*this)[1][1] -
(*this)[0][1] * (*this)[1][0];
}
// [ a, b, c ]
// Given matrix A = [ d, e, f ]
// [ g, h, i ]
//
// det (A) = aei + bfg + cdh - afg - bdi - ceg
// det (A) = a(ei - fg) + b(fg - di) + c(dh - eg)
double
matrix::determinant3x3 (void) const {
CHECK_EQ (m_rows, 3);
CHECK_EQ (m_columns, 3);
return (*this)[0][0] * (*this)[1][1] * (*this)[2][2] + // aei
(*this)[0][1] * (*this)[1][2] * (*this)[2][0] + // bfg
(*this)[0][2] * (*this)[1][0] * (*this)[2][1] - // cdh
(*this)[0][0] * (*this)[1][2] * (*this)[2][1] - // afh
(*this)[0][1] * (*this)[1][0] * (*this)[2][2] - // bdi
(*this)[0][2] * (*this)[1][1] * (*this)[2][0]; // ceg
}
// From libMathematics, http://www.geometrictools.com/
double
matrix::determinant4x4 (void) const {
CHECK_EQ (m_rows, 4);
CHECK_EQ (m_columns, 4);
double a0 = m_data[ 0] * m_data[ 5] - m_data[ 1] * m_data[ 4],
a1 = m_data[ 0] * m_data[ 6] - m_data[ 2] * m_data[ 4],
a2 = m_data[ 0] * m_data[ 7] - m_data[ 3] * m_data[ 4],
a3 = m_data[ 1] * m_data[ 6] - m_data[ 2] * m_data[ 5],
a4 = m_data[ 1] * m_data[ 7] - m_data[ 3] * m_data[ 5],
a5 = m_data[ 2] * m_data[ 7] - m_data[ 3] * m_data[ 6],
b0 = m_data[ 8] * m_data[13] - m_data[ 9] * m_data[12],
b1 = m_data[ 8] * m_data[14] - m_data[10] * m_data[12],
b2 = m_data[ 8] * m_data[15] - m_data[11] * m_data[12],
b3 = m_data[ 9] * m_data[14] - m_data[10] * m_data[13],
b4 = m_data[ 9] * m_data[15] - m_data[11] * m_data[13],
b5 = m_data[10] * m_data[15] - m_data[11] * m_data[14];
return a0 * b5 - a1 * b4 + a2 * b3 + a3 * b2 - a4 * b1 + a5 * b0;
}
matrix
matrix::inverse (void) const {
if (m_rows != m_columns)
not_implemented ();
switch (m_rows) {
case 2: return inverse2x2 ();
case 3: return inverse3x3 ();
case 4: return inverse4x4 ();
}
not_implemented ();
}
matrix
matrix::inverse2x2 (void) const {
CHECK (m_rows == 2);
CHECK (m_columns == 2);
double det = determinant2x2 ();
if (almost_equal (det, 0.))
throw not_invertible ();
return matrix (2, 2, { (*this)[1][1], -(*this)[0][1],
-(*this)[1][0], (*this)[0][0] }) /= det;
}
// [ a, b, c ]
// Given matrix A = [ d, e, f ]
// [ g, h, i ]
//
matrix
matrix::inverse3x3 (void) const {
CHECK (m_rows == 3);
CHECK (m_columns == 3);
double det = determinant3x3();
if (almost_equal (det, 0.))
throw not_invertible ();
matrix val (m_rows, m_columns, {
(*this)[1][1] * (*this)[2][2] - (*this)[1][2] * (*this)[2][1], // ei - fh
(*this)[0][2] * (*this)[2][1] - (*this)[0][1] * (*this)[2][2], // ch - bi
(*this)[0][1] * (*this)[1][2] - (*this)[0][2] * (*this)[1][1], // bf - ce
(*this)[1][2] * (*this)[2][0] - (*this)[1][0] * (*this)[2][2], // fg - di
(*this)[0][0] * (*this)[2][2] - (*this)[0][2] * (*this)[2][0], // ai - cg
(*this)[0][2] * (*this)[1][0] - (*this)[0][0] * (*this)[1][2], // cd - af
(*this)[1][0] * (*this)[2][1] - (*this)[1][1] * (*this)[2][0], // dh - eg
(*this)[0][1] * (*this)[2][0] - (*this)[0][0] * (*this)[2][1], // bg - ah
(*this)[0][0] * (*this)[1][1] - (*this)[0][1] * (*this)[1][0] // ae - bd
});
return val /= det;
//matrix val ({ vector::cross ((*this)[1], (*this)[2], 3),
// vector::cross ((*this)[2], (*this)[0], 3),
// vector::cross ((*this)[0], (*this)[1], 3) });
//return val /= determinant3x3 ();
}
matrix
matrix::inverse4x4 (void) const {
double a0 = m_data[ 0] * m_data[ 5] - m_data[ 1] * m_data[ 4],
a1 = m_data[ 0] * m_data[ 6] - m_data[ 2] * m_data[ 4],
a2 = m_data[ 0] * m_data[ 7] - m_data[ 3] * m_data[ 4],
a3 = m_data[ 1] * m_data[ 6] - m_data[ 2] * m_data[ 5],
a4 = m_data[ 1] * m_data[ 7] - m_data[ 3] * m_data[ 5],
a5 = m_data[ 2] * m_data[ 7] - m_data[ 3] * m_data[ 6],
b0 = m_data[ 8] * m_data[13] - m_data[ 9] * m_data[12],
b1 = m_data[ 8] * m_data[14] - m_data[10] * m_data[12],
b2 = m_data[ 8] * m_data[15] - m_data[11] * m_data[12],
b3 = m_data[ 9] * m_data[14] - m_data[10] * m_data[13],
b4 = m_data[ 9] * m_data[15] - m_data[11] * m_data[13],
b5 = m_data[10] * m_data[15] - m_data[11] * m_data[14];
double det = a0 * b5 - a1 * b4 + a2 * b3 + a3 * b2 - a4 * b1 + a5 * b0;
if (almost_equal (det, 0.))
throw not_invertible ();
return matrix (4, 4, {
+ m_data[ 5] * b5 - m_data[ 6] * b4 + m_data[ 7] * b3,
- m_data[ 1] * b5 + m_data[ 2] * b4 - m_data[ 3] * b3,
+ m_data[13] * a5 - m_data[14] * a4 + m_data[15] * a3,
- m_data[ 9] * a5 + m_data[10] * a4 - m_data[11] * a3,
- m_data[ 4] * b5 + m_data[ 6] * b2 - m_data[ 7] * b1,
+ m_data[ 0] * b5 - m_data[ 2] * b2 + m_data[ 3] * b1,
- m_data[12] * a5 + m_data[14] * a2 - m_data[15] * a1,
+ m_data[ 8] * a5 - m_data[10] * a2 + m_data[11] * a1,
+ m_data[ 4] * b4 - m_data[ 5] * b2 + m_data[ 7] * b0,
- m_data[ 0] * b4 + m_data[ 1] * b2 - m_data[ 3] * b0,
+ m_data[12] * a4 - m_data[13] * a2 + m_data[15] * a0,
- m_data[ 8] * a4 + m_data[ 9] * a2 - m_data[11] * a0,
- m_data[ 4] * b3 + m_data[ 5] * b1 - m_data[ 6] * b0,
+ m_data[ 0] * b3 - m_data[ 1] * b1 + m_data[ 2] * b0,
- m_data[12] * a3 + m_data[13] * a1 - m_data[14] * a0,
+ m_data[ 8] * a3 - m_data[ 9] * a1 + m_data[10] * a0
}) /= det;
}
matrix
matrix::zeroes (size_t diag)
{ return zeroes (diag, diag); }
matrix
matrix::zeroes (size_t rows, size_t columns) {
matrix m (rows, columns);
std::fill (m.m_data.get (), m.m_data.get () + m.size (), 0.0);
return m;
}
matrix
matrix::identity (size_t diag) {
matrix val (zeroes (diag));
for (unsigned int i = 0; i < diag; ++i)
val[i][i] = 1.0;
return val;
}
matrix
matrix::magic (size_t n) {
CHECK_GT (n, 2);
if (n % 2 == 1)
return magic_odd (n);
if (n % 4 == 0)
return magic_even_single (n);
return magic_even_double (n);
}
// Use the 'siamese' method. Start from the top centre, progress up-left one.
// If filled then drop down one row instead. Wrap around indexing.
matrix
matrix::magic_odd (size_t n) {
CHECK_GT (n, 2);
CHECK_EQ (n % 2, 1);
matrix val (zeroes (n));
for (unsigned int i = 1, x = n / 2, y = 0; i <= n * n; ++i) {
val[y][x] = i;
unsigned int x1 = (x + 1) % n,
y1 = (y + n - 1) % n;
if (!almost_equal (val[y1][x1], 0)) {
x1 = x;
y1 = (y + 1) % n;
}
x = x1;
y = y1;
}
return val;
}
matrix
matrix::magic_even_single (size_t)
{ not_implemented (); }
matrix
matrix::magic_even_double (size_t)
{ not_implemented (); }

117
maths/matrix.hpp
View File

@ -1,117 +0,0 @@
/*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
* Copyright 2010 Danny Robson <danny@nerdcruft.net>
*/
#ifndef __UTIL_MATHS_MATRIX_HPP
#define __UTIL_MATHS_MATRIX_HPP
#include "vector.hpp"
#include <algorithm>
#include <assert.h>
#include <initializer_list>
#include <iostream>
#include <memory>
#include <stdexcept>
namespace maths {
class matrix {
protected:
size_t m_rows,
m_columns;
std::unique_ptr<double[]> m_data;
public:
matrix (size_t _rows, size_t _columns);
matrix (size_t _rows,
size_t _columns,
const std::initializer_list <double> &_data);
matrix (const std::initializer_list <vector> &_data);
matrix (const matrix &rhs);
matrix (matrix &&rhs);
~matrix();
void sanity (void) const;
const double * operator [] (unsigned int row) const;
double * operator [] (unsigned int row);
const double * data (void) const;
matrix& operator =(const matrix &rhs);
matrix operator * (double scalar) const;
matrix& operator *=(double scalar);
matrix operator * (const matrix &rhs) const;
matrix& operator *=(const matrix &rhs);
matrix& operator /=(double scalar);
matrix operator + (double scalar) const;
matrix& operator +=(double scalar);
matrix& operator -=(double scalar);
bool operator ==(const matrix &rhs) const;
//matrix transpose (void) const { ; }
size_t rows (void) const;
size_t columns (void) const;
size_t size (void) const;
/// Checks if this is a sqaure matrix, with a zero final column
/// and row (excepting the final diagonal entry).
bool is_homogeneous (void) const;
bool is_square (void) const;
bool is_magic (void) const;
public:
double determinant (void) const;
matrix inverse (void) const;
protected:
double determinant2x2 (void) const;
double determinant3x3 (void) const;
double determinant4x4 (void) const;
matrix inverse2x2 (void) const;
matrix inverse3x3 (void) const;
matrix inverse4x4 (void) const;
public:
static matrix zeroes (size_t n);
static matrix zeroes (size_t rows, size_t columns);
static matrix identity (size_t n);
/// Generate a magic square of order 'n'
static matrix magic (size_t n);
protected:
/// Generate a magic square with 'n' odd
static matrix magic_odd (size_t n);
/// Generate a magic square with 'n' divisible by 2, and not 4
static matrix magic_even_single (size_t n);
/// Generate a magic square with 'n' divisible by 4, and not 2
static matrix magic_even_double (size_t n);
};
class not_invertible : public std::runtime_error {
public:
not_invertible ():
std::runtime_error ("not_invertible")
{ ; }
};
}
#endif

81
maths/vector.cpp
View File

@ -1,81 +0,0 @@
#include "vector.hpp"
#include "debug.hpp"
#include <numeric>
using namespace maths;
/* Constructors */
vector::vector (const std::initializer_list<double> &_data):
m_data (_data)
{ ; }
vector::vector (unsigned int _size):
m_data (_size)
{ ; }
vector::vector (const double *restrict _data,
unsigned int _size):
m_data (_size)
{ std::copy (_data, _data + _size, m_data.begin ()); }
vector::vector (const vector &rhs):
m_data (rhs.m_data)
{ ; }
vector::vector (const vector &&rhs):
m_data (std::move (rhs.m_data))
{ ; }
vector::~vector (void)
{ ; }
/* element accessors */
const double*
vector::data (void) const
{ return &m_data[0]; }
double &
vector::operator[] (unsigned int offset)
{ return m_data[offset]; }
const double&
vector::operator[] (unsigned int offset) const
{ return m_data[offset]; }
unsigned int
vector::size (void) const
{ return m_data.size (); }
/* dot and cross products */
double vector::dot (const double *restrict A,
const double *restrict B,
unsigned int size)
{ return std::inner_product(A, A + size, B, 0.0); }
vector vector::cross (const double *restrict A,
const double *restrict B,
unsigned int size) {
CHECK_EQ (size, 3);
(void)size;
return vector ({ A[1] * B[2] - A[2] * B[1],
A[2] * B[0] - A[0] * B[2],
A[0] * B[1] - A[1] * B[0] });
}

58
maths/vector.hpp
View File

@ -1,58 +0,0 @@
/*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
* Copyright 2011 Danny Robson <danny@nerdcruft.net>
*/
#ifndef __UTIL_MATHS_VECTOR_HPP
#define __UTIL_MATHS_VECTOR_HPP
#include <vector>
#include <initializer_list>
namespace maths {
class vector {
protected:
std::vector<double> m_data;
public:
vector (const std::initializer_list<double> &_data);
explicit
vector (unsigned int _size);
vector (const double *restrict _data,
unsigned int _size);
vector (const vector &rhs);
vector (const vector &&rhs);
~vector (void);
const double* data (void) const;
double& operator[] (unsigned int);
const double& operator[] (unsigned int) const;
unsigned int size (void) const;
static double dot (const double *restrict A,
const double *restrict B,
unsigned int size);
static vector cross (const double *restrict A,
const double *restrict B,
unsigned int size);
};
}
#endif

136
test/maths_matrix.cpp
View File

@ -1,136 +0,0 @@
#include "maths/matrix.hpp"
#include "debug.hpp"
#include "maths.hpp"
#include "tap.hpp"
#include <iostream>
#include <cmath>
using namespace maths;
using namespace std;
std::ostream&
operator <<(std::ostream &os, const matrix &m) {
for (unsigned int i = 0; i < m.rows (); ++i) {
for (unsigned int j = 0; j < m.columns (); ++j) {
os << m[i][j];
if (j != m.columns () - 1)
os << ", ";
}
if (i != m.rows () - 1)
os << "\n";
}
return os;
}
void
test_zeroes (const matrix &m) {
assert (m.rows ());
assert (m.columns ());
for (unsigned int i = 0; i < m.rows (); ++i)
for (unsigned int j = 0; j < m.columns (); ++j)
CHECK (almost_equal (m[i][j], 0.0));
}
void
test_identity (const matrix &m) {
assert (m.rows () == m.columns ());
for (unsigned int i = 0; i < m.rows (); ++i)
for (unsigned int j = 0; j < m.columns (); ++j)
if (i == j) {
CHECK (almost_equal (m[i][j], 1.0));
} else {
CHECK (almost_equal (m[i][j], 0.0));
}
}
int
main (void) {
for (unsigned int i = 1; i < 10; ++i) {
test_zeroes (matrix::zeroes (i));
test_identity (matrix::identity (i));
}
for (unsigned int i = 3; i < 10; i += 2)
CHECK (matrix::magic (i).is_magic ());
// Create a small matrix with unique element values for comparison tests.
// This should be non-square so that row- vs. column-major problems can
// be seen.
matrix a4x2 (4, 2, { 0, 1,
2, 3,
4, 5,
6, 7 });
CHECK_EQ (a4x2, a4x2);
// Test that copy constructors work correctly. Keep this value around so
// that we can check the following operators don't modify the original
// value.
CHECK_EQ (a4x2, a4x2);
// Check multiplication by identity results in the original value.
CHECK_EQ (a4x2, a4x2 * matrix::identity (a4x2.columns ()));
matrix seq2x2(2, 2, { 1, 2, 3, 4 });
matrix magic3(3, 3, { 2, 7, 6,
9, 5, 1,
4, 3, 8 });
matrix magic4(4, 4, { 16, 2, 3, 13,
5, 11, 10, 8,
9, 7, 6, 12,
4, 14, 15, 1 });
CHECK_EQ (magic3[0][0], 2.0);
CHECK_EQ (magic3[0][1], 7.0);
CHECK_EQ (magic3[0][2], 6.0);
CHECK_EQ (magic3[1][0], 9.0);
CHECK_EQ (magic3[1][1], 5.0);
CHECK_EQ (magic3[1][2], 1.0);
CHECK_EQ (magic3[2][0], 4.0);
CHECK_EQ (magic3[2][1], 3.0);
CHECK_EQ (magic3[2][2], 8.0);
CHECK_EQ (seq2x2.determinant (), -2.0);
CHECK_EQ (magic3.determinant (), -360.0);
CHECK ( seq2x2.is_square ());
CHECK ( magic3.is_square ());
CHECK (! a4x2.is_square ());
CHECK_EQ (seq2x2.inverse (), matrix (2, 2, { -2.0, 1.0,
1.5, -0.5 }));
CHECK_EQ (magic3.inverse (), matrix (3, 3, { -37.0, 38.0, 23.0,
68.0, 8.0, -52.0,
- 7.0, -22.0, 53.0 }) /= 360.0);
matrix invertible4 (4, 4, { 4, 14, 15, 1,
9, 7, 6, 12,
5, 11, 10, 8,
0, 0, 0, 1 });
CHECK_EQ (invertible4.inverse (), matrix (4, 4, { 4, 25, -21, -136,
-60, -35, 111, -408,
64, 26, -98, 408,
0, 0, 0, 136 }) /= 136);
const matrix homo3x3 (3, 3, { 1, 2, 0,
3, 4, 0,
0, 0, 1 });
CHECK (homo3x3.is_homogeneous ());
CHECK (!matrix::zeroes (3).is_homogeneous ());
CHECK ( matrix::identity (3).is_homogeneous ());
CHECK (invertible4.is_homogeneous ());
util::TAP::logger tap;
tap.skip ("convert to TAP");
}