Big math::vector,point,matrix refactor: fully support dynamic (runtime) size in basic operations

This commit is contained in:
Nikita Lisitsa 2025-12-04 17:22:32 +03:00
parent a8e08bf6d3
commit a5acb9534b
6 changed files with 442 additions and 232 deletions

View file

@ -2,7 +2,6 @@
#include <psemek/math/dynamic.hpp>
#include <psemek/util/exception.hpp>
#include <psemek/util/array.hpp>
namespace psemek::math::detail
{
@ -15,10 +14,42 @@ namespace psemek::math::detail
{}
};
template <std::size_t N>
struct dynamic_dimensions
{};
template <>
struct dynamic_dimensions<dynamic>
{
std::size_t size = 0;
};
template <std::size_t N>
void check_dynamic_size(dynamic_dimensions<N> d1, dynamic_dimensions<N> d2)
{
if constexpr (N == dynamic)
{
if (d1.size != d2.size)
throw dynamic_size_mismatch(d1.size, d2.size);
}
}
template <typename T, std::size_t N>
struct array
{
using type = T[N];
struct type
{
static constexpr std::size_t size = N;
T data[N];
type(dynamic_dimensions<N>){}
dynamic_dimensions<N> dimensions() const { return {}; }
T const & operator[](std::size_t i) const { return data[i]; }
T & operator[](std::size_t i) { return data[i]; }
};
};
template <typename T>
@ -26,17 +57,35 @@ namespace psemek::math::detail
{
struct type
{
static constexpr std::size_t size = 0;
type(dynamic_dimensions<0>){}
dynamic_dimensions<0> dimensions() const { return {}; }
T const & operator[](std::size_t) const { throw empty_array_exception{}; }
T & operator[](std::size_t) { throw empty_array_exception{}; }
// type operator + (std::size_t) const { return *this; }
};
};
template <typename T>
struct array<T, dynamic>
{
using type = util::array<T>;
struct type
{
std::size_t size;
std::unique_ptr<T[]> data;
type(dynamic_dimensions<dynamic> dimensions)
: size(dimensions.size)
, data(std::make_unique_for_overwrite<T[]>(dimensions.size))
{}
dynamic_dimensions<dynamic> dimensions() const { return {.size = size}; }
T const & operator[](std::size_t i) const { return data[i]; }
T & operator[](std::size_t i) { return data[i]; }
};
};
template <typename T, typename ... Args>

View file

@ -5,10 +5,37 @@
namespace psemek::math::detail
{
template <std::size_t R, std::size_t C>
struct dynamic_dimensions_2d
{
dynamic_dimensions<R> rows = {};
dynamic_dimensions<C> columns = {};
};
template <std::size_t R, std::size_t C>
void check_dynamic_size(dynamic_dimensions_2d<R, C> d1, dynamic_dimensions_2d<R, C> d2)
{
check_dynamic_size(d1.rows, d2.rows);
check_dynamic_size(d1.columns, d2.columns);
}
template <typename T, std::size_t R, std::size_t C>
struct array_2d
{
using type = T[R][C];
struct type
{
T data[R * C];
static constexpr std::size_t rows = R;
static constexpr std::size_t columns = C;
type(dynamic_dimensions_2d<R, C>) {}
dynamic_dimensions_2d<R, C> dimensions() const { return {}; }
T const * operator[](std::size_t row) const { return data + C * row; }
T * operator[](std::size_t row) { return data + C * row; }
};
};
template <typename T, std::size_t R>
@ -19,6 +46,10 @@ namespace psemek::math::detail
static constexpr std::size_t rows = R;
static constexpr std::size_t columns = 0;
type(dynamic_dimensions_2d<R, 0>) {}
dynamic_dimensions_2d<R, 0> dimensions() const { return {}; }
T const * operator[](std::size_t) const { throw empty_array_exception{}; }
T * operator[](std::size_t) { throw empty_array_exception{}; }
};
@ -32,10 +63,12 @@ namespace psemek::math::detail
std::size_t rows;
static constexpr std::size_t columns = 0;
type(std::size_t rows)
: rows(rows)
type(dynamic_dimensions_2d<dynamic, 0> dimensions)
: rows(dimensions.rows.size)
{}
dynamic_dimensions_2d<dynamic, 0> dimensions() const { return {.rows = rows}; }
T const * operator[](std::size_t) const { throw empty_array_exception{}; }
T * operator[](std::size_t) { throw empty_array_exception{}; }
};
@ -49,6 +82,10 @@ namespace psemek::math::detail
static constexpr std::size_t rows = 0;
static constexpr std::size_t columns = C;
type(dynamic_dimensions_2d<0, C>) {}
dynamic_dimensions_2d<0, C> dimensions() const { return {}; }
T const * operator[](std::size_t) const { throw empty_array_exception{}; }
T * operator[](std::size_t) { throw empty_array_exception{}; }
};
@ -62,10 +99,12 @@ namespace psemek::math::detail
static constexpr std::size_t rows = 0;
std::size_t columns;
type(std::size_t columns)
: columns(columns)
type(dynamic_dimensions_2d<0, dynamic> dimensions)
: columns(dimensions.columns.size)
{}
dynamic_dimensions_2d<0, dynamic> dimensions() const { return {.columns = columns}; }
T const * operator[](std::size_t) const { throw empty_array_exception{}; }
T * operator[](std::size_t) { throw empty_array_exception{}; }
};
@ -79,6 +118,10 @@ namespace psemek::math::detail
static constexpr std::size_t rows = 0;
static constexpr std::size_t columns = 0;
type(dynamic_dimensions_2d<0, 0>) {}
dynamic_dimensions_2d<0, 0> dimensions() const { return {}; }
T const * operator[](std::size_t) const { throw empty_array_exception{}; }
T * operator[](std::size_t) { throw empty_array_exception{}; }
};
@ -93,20 +136,12 @@ namespace psemek::math::detail
std::size_t columns;
std::unique_ptr<T[]> data;
type()
: columns(0)
{}
type(std::size_t columns)
: columns(columns)
type(dynamic_dimensions_2d<R, dynamic> dimensions)
: columns(dimensions.columns.size)
, data(std::make_unique_for_overwrite<T[]>(rows * columns))
{}
type([[maybe_unused]] std::size_t rows, std::size_t columns)
: type(columns)
{
assert(rows == R);
}
dynamic_dimensions_2d<R, dynamic> dimensions() const { return {.columns = columns}; }
T * operator[] (std::size_t row)
{
@ -138,20 +173,12 @@ namespace psemek::math::detail
static constexpr std::size_t columns = C;
std::unique_ptr<T[]> data;
type()
: rows(0)
{}
type(std::size_t rows)
: rows(rows)
type(dynamic_dimensions_2d<dynamic, C> dimensions)
: rows(dimensions.rows.size)
, data(std::make_unique_for_overwrite<T[]>(rows * columns))
{}
type(std::size_t rows, [[maybe_unused]] std::size_t columns)
: type(rows)
{
assert(columns == C);
}
dynamic_dimensions_2d<dynamic, C> dimensions() const { return {.rows = rows}; }
T * operator[] (std::size_t row)
{
@ -183,17 +210,14 @@ namespace psemek::math::detail
std::size_t columns;
std::unique_ptr<T[]> data;
type()
: rows(0)
, columns(0)
{}
type(std::size_t rows, std::size_t columns)
: rows(rows)
, columns(columns)
type(dynamic_dimensions_2d<dynamic, dynamic> dimensions)
: rows(dimensions.rows.size)
, columns(dimensions.columns.size)
, data(std::make_unique_for_overwrite<T[]>(rows * columns))
{}
dynamic_dimensions_2d<dynamic, dynamic> dimensions() const { return {.rows = rows, .columns = columns}; }
T * operator[] (std::size_t row)
{
return data.get() + row * columns;

View file

@ -18,10 +18,4 @@ namespace psemek::math
{}
};
inline void check_dynamic_size(std::size_t size1, std::size_t size2)
{
if (size1 != size2)
throw dynamic_size_mismatch(size1, size2);
}
}

View file

@ -19,31 +19,31 @@ namespace psemek::math
static constexpr std::size_t static_rows = R;
static constexpr std::size_t static_columns = C;
using dynamic_dimensions_type = detail::dynamic_dimensions_2d<R, C>;
typename detail::array_2d<T, R, C>::type coords;
// Internal
matrix(dynamic_dimensions_type dimensions);
matrix();
template <typename ... Args>
matrix(Args && ... args) requires (R != dynamic && C != dynamic && sizeof...(Args) == R * C && detail::all_convertible_to<T, Args...>::value);
explicit matrix(std::size_t rows) requires(R == dynamic && C != dynamic);
explicit matrix(std::size_t columns) requires(R != dynamic && C == dynamic);
matrix(std::size_t rows, std::size_t columns) requires(R == dynamic && C == dynamic);
std::size_t rows() const
{
if constexpr (R == dynamic)
{
return coords.rows;
}
else
{
return R;
}
}
std::size_t columns() const
{
if constexpr (C == dynamic)
{
return coords.columns;
}
else
{
return C;
}
}
auto operator[](std::size_t i)
{
@ -85,32 +85,72 @@ namespace psemek::math
matrix & operator += (matrix const & v);
matrix & operator -= (matrix const & v);
static matrix zero();
static matrix identity();
static matrix scalar(T const & s);
static matrix zero(dynamic_dimensions_type dimensions = {});
static matrix identity(dynamic_dimensions_type dimensions = {});
static matrix scalar(T const & s, dynamic_dimensions_type dimensions = {});
static matrix diagonal(vector<T, std::min(R, C)> const & d);
};
template <typename T, std::size_t R, std::size_t C>
matrix<T, R, C> matrix<T, R, C>::zero()
void check_dynamic_size(matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
{
matrix<T, R, C> m;
check_dynamic_size(m1.coords.dimensions(), m2.coords.dimensions());
}
template <typename T, std::size_t R, std::size_t C>
matrix<T, R, C>::matrix(dynamic_dimensions_type dimensions)
: coords(dimensions)
{}
template <typename T, std::size_t R, std::size_t C>
matrix<T, R, C>::matrix()
: matrix(dynamic_dimensions_type{})
{}
template <typename T, std::size_t R, std::size_t C>
template <typename ... Args>
matrix<T, R, C>::matrix(Args && ... args) requires (R != dynamic && C != dynamic && sizeof...(Args) == R * C && detail::all_convertible_to<T, Args...>::value)
: coords({})
{
auto out = values().begin();
((*out++ = args), ...);
}
template <typename T, std::size_t R, std::size_t C>
matrix<T, R, C>::matrix(std::size_t rows) requires(R == dynamic && C != dynamic)
: coords({.rows = rows})
{}
template <typename T, std::size_t R, std::size_t C>
matrix<T, R, C>::matrix(std::size_t columns) requires(R != dynamic && C == dynamic)
: coords({.columns = columns})
{}
template <typename T, std::size_t R, std::size_t C>
matrix<T, R, C>::matrix(std::size_t rows, std::size_t columns) requires(R == dynamic && C == dynamic)
: coords({.rows = rows, .columns = columns})
{}
template <typename T, std::size_t R, std::size_t C>
matrix<T, R, C> matrix<T, R, C>::zero(dynamic_dimensions_type dimensions)
{
matrix<T, R, C> m(dimensions);
for (auto & v : m.values())
v = T(0);
return m;
}
template <typename T, std::size_t R, std::size_t C>
matrix<T, R, C> matrix<T, R, C>::identity()
matrix<T, R, C> matrix<T, R, C>::identity(dynamic_dimensions_type dimensions)
{
return scalar(T(1));
return scalar(T(1), dimensions);
}
template <typename T, std::size_t R, std::size_t C>
matrix<T, R, C> matrix<T, R, C>::scalar(T const & s)
matrix<T, R, C> matrix<T, R, C>::scalar(T const & s, dynamic_dimensions_type dimensions)
{
matrix<T, R, C> m = zero();
for (std::size_t i = 0; i < std::min(R, C); ++i)
matrix<T, R, C> m = zero(dimensions);
for (std::size_t i = 0; i < std::min(m.rows(), m.columns()); ++i)
m[i][i] = s;
return m;
}
@ -118,7 +158,17 @@ namespace psemek::math
template <typename T, std::size_t R, std::size_t C>
matrix<T, R, C> matrix<T, R, C>::diagonal(vector<T, std::min(R, C)> const & d)
{
matrix<T, R, C> m = zero();
dynamic_dimensions_type dimensions;
if constexpr (R == dynamic)
{
dimensions.rows = d.dimension();
}
if constexpr (C == dynamic)
{
dimensions.columns = d.dimension();
}
matrix<T, R, C> m = zero(dimensions);
for (std::size_t i = 0; i < std::min(R, C); ++i)
m[i][i] = d[i];
return m;
@ -127,25 +177,33 @@ namespace psemek::math
template <typename T1, typename T, std::size_t R, std::size_t C>
matrix<T1, R, C> cast(matrix<T, R, C> const & m)
{
matrix<T1, R, C> r;
for (std::size_t i = 0; i < R; ++i)
for (std::size_t j = 0; j < C; ++j)
r[i][j] = static_cast<T1>(m[i][j]);
matrix<T1, R, C> r(m.coords.dimensions());
auto out = r.values().begin();
for (auto const & value : m.values())
*out++ = static_cast<T1>(value);
return r;
}
template <typename T, std::size_t R, std::size_t C>
std::strong_ordering operator <=> (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
{
for (std::size_t i = 0; i < R; ++i)
check_dynamic_size(m1, m2);
auto begin1 = m1.values().begin();
auto begin2 = m2.values().begin();
auto end1 = m1.values().end();
// Intentionally don't use std::lexicographical_compare_three_way to fake strong ordering
// (actual ordering might be partial & weak due to NaN's and signed zeros
while (begin1 != end1)
{
for (std::size_t j = 0; j < C; ++j)
{
if (m1[i][j] < m2[i][j])
if (*begin1 < *begin2)
return std::strong_ordering::less;
else if (m1[i][j] > m2[i][j])
else if (*begin1 > *begin2)
return std::strong_ordering::greater;
}
++begin1;
++begin2;
}
return std::strong_ordering::equal;
}
@ -189,30 +247,30 @@ namespace psemek::math
template <typename T, std::size_t R, std::size_t C>
matrix<T, R, C> operator * (matrix<T, R, C> const & m, T const & s)
{
matrix<T, R, C> r;
for (std::size_t i = 0; i < R; ++i)
for (std::size_t j = 0; j < C; ++j)
r[i][j] = m[i][j] * s;
matrix<T, R, C> r(m.coords.dimensions());
auto out = r.values().begin();
for (auto const & value : m.values())
*out++ = value * s;
return r;
}
template <typename T, std::size_t R, std::size_t C>
matrix<T, R, C> operator * (T const & s, matrix<T, R, C> const & m)
{
matrix<T, R, C> r;
for (std::size_t i = 0; i < R; ++i)
for (std::size_t j = 0; j < C; ++j)
r[i][j] = m[i][j] * s;
matrix<T, R, C> r(m.coords.dimensions());
auto out = r.values().begin();
for (auto const & value : m.values())
*out++ = s * value;
return r;
}
template <typename T, std::size_t R, std::size_t C>
matrix<T, R, C> operator / (matrix<T, R, C> const & m, T const & s)
{
matrix<T, R, C> r;
for (std::size_t i = 0; i < R; ++i)
for (std::size_t j = 0; j < C; ++j)
r[i][j] = m[i][j] / s;
matrix<T, R, C> r(m.coords.dimensions());
auto out = r.values().begin();
for (auto const & value : m.values())
*out++ = value / s;
return r;
}
@ -233,30 +291,34 @@ namespace psemek::math
template <typename T, std::size_t R, std::size_t C>
matrix<T, R, C> operator - (matrix<T, R, C> const & m)
{
matrix<T, R, C> r;
for (std::size_t i = 0; i < R; ++i)
for (std::size_t j = 0; j < C; ++j)
r[i][j] = -m[i][j];
matrix<T, R, C> r(m.coords.dimensions());
auto out = r.values().begin();
for (auto const & value : m.values())
*out++ = -value;
return r;
}
template <typename T, std::size_t R, std::size_t C>
matrix<T, R, C> operator + (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
{
matrix<T, R, C> r;
for (std::size_t i = 0; i < R; ++i)
for (std::size_t j = 0; j < C; ++j)
r[i][j] = m1[i][j] + m2[i][j];
check_dynamic_size(m1, m2);
matrix<T, R, C> r(m1.coords.dimensions());
auto in1 = m1.values().begin();
auto in2 = m2.values().begin();
for (auto & value : r.values())
value = (*in1++) + (*in2++);
return r;
}
template <typename T, std::size_t R, std::size_t C>
matrix<T, R, C> operator - (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
{
matrix<T, R, C> r;
for (std::size_t i = 0; i < R; ++i)
for (std::size_t j = 0; j < C; ++j)
r[i][j] = m1[i][j] - m2[i][j];
check_dynamic_size(m1, m2);
matrix<T, R, C> r(m1.coords.dimensions());
auto in1 = m1.values().begin();
auto in2 = m2.values().begin();
for (auto & value : r.values())
value = (*in1++) - (*in2++);
return r;
}
@ -277,11 +339,11 @@ namespace psemek::math
template <typename T, std::size_t R, std::size_t C>
vector<T, R> operator * (matrix<T, R, C> const & m, vector<T, C> const & v)
{
vector<T, R> r;
for (std::size_t i = 0; i < R; ++i)
detail::check_dynamic_size(m.coords.dimensions().columns, v.coords.dimensions());
auto r = vector<T, R>::zero(m.coords.dimensions().rows);
for (std::size_t i = 0; i < m.rows(); ++i)
{
r[i] = T{};
for (std::size_t j = 0; j < C; ++j)
for (std::size_t j = 0; j < m.columns(); ++j)
r[i] += m[i][j] * v[j];
}
return r;
@ -290,12 +352,12 @@ namespace psemek::math
template <typename T, std::size_t R, std::size_t C>
vector<T, C> operator * (vector<T, R> const & v, matrix<T, R, C> const & m)
{
vector<T, C> r;
for (std::size_t j = 0; j < C; ++j)
detail::check_dynamic_size(v.coords.dimensions(), m.coords.dimensions().rows);
auto r = vector<T, C>::zero(m.coords.dimensions().columns);
for (std::size_t i = 0; i < m.rows(); ++i)
{
r[j] = T{};
for (std::size_t i = 0; i < R; ++i)
r[i] += v[j] * m[i][j];
for (std::size_t j = 0; j < m.columns(); ++j)
r[j] += v[i] * m[i][j];
}
return r;
}
@ -303,14 +365,13 @@ namespace psemek::math
template <typename T, std::size_t R, std::size_t K, std::size_t C>
matrix<T, R, C> operator * (matrix<T, R, K> const & m1, matrix<T, K, C> const & m2)
{
matrix<T, R, C> r;
for (std::size_t i = 0; i < R; ++i)
detail::check_dynamic_size(m1.coords.dimensions().columns, m2.coords.dimensions().rows);
auto r = matrix<T, R, C>::zero();
for (std::size_t i = 0; i < m1.rows(); ++i)
{
for (std::size_t j = 0; j < C; ++j)
for (std::size_t j = 0; j < m2.columns(); ++j)
{
r[i][j] = T{};
for (std::size_t k = 0; k < K; ++k)
for (std::size_t k = 0; k < m1.columns(); ++k)
r[i][j] += m1[i][k] * m2[k][j];
}
}
@ -320,9 +381,9 @@ namespace psemek::math
template <typename T, std::size_t R, std::size_t C>
matrix<T, C, R> transpose(matrix<T, R, C> const & m)
{
matrix<T, C, R> r;
for (std::size_t i = 0; i < R; ++i)
for (std::size_t j = 0; j < C; ++j)
matrix<T, C, R> r({.rows = m.coords.dimensions().columns, .columns = m.coords.dimensions().rows});
for (std::size_t i = 0; i < m.rows(); ++i)
for (std::size_t j = 0; j < m.columns(); ++j)
r[j][i] = m[i][j];
return r;
}
@ -330,9 +391,11 @@ namespace psemek::math
template <typename T, std::size_t N, typename ... Rows>
auto by_rows(vector<T, N> const & v, Rows const & ... rows)
{
(check_dynamic_size(v, rows), ...);
vector<T, N> m[] = {v, rows...};
matrix<T, sizeof...(Rows) + 1, N> r;
matrix<T, sizeof...(Rows) + 1, N> r({.columns = v.coords.dimensions()});
for (std::size_t i = 0; i < r.rows(); ++i)
for (std::size_t j = 0; j < r.columns(); ++j)
r[i][j] = m[i][j];
@ -342,9 +405,10 @@ namespace psemek::math
template <typename T, std::size_t N, typename ... Columns>
auto by_columns(vector<T, N> const & v, Columns const & ... columns)
{
(check_dynamic_size(v, columns), ...);
vector<T, N> m[] = {v, columns...};
matrix<T, N, sizeof...(Columns) + 1> r;
matrix<T, N, sizeof...(Columns) + 1> r({.rows = v.coords.dimensions()});
for (std::size_t i = 0; i < r.rows(); ++i)
for (std::size_t j = 0; j < r.columns(); ++j)
r[i][j] = m[j][i];
@ -354,8 +418,8 @@ namespace psemek::math
template <typename T, std::size_t R, std::size_t C>
vector<T, C> row(matrix<T, R, C> const & m, std::size_t i)
{
vector<T, C> r;
for (std::size_t j = 0; j < C; ++j)
vector<T, C> r(m.coords.dimensions().columns);
for (std::size_t j = 0; j < m.columns(); ++j)
r[j] = m[i][j];
return r;
}
@ -363,8 +427,8 @@ namespace psemek::math
template <typename T, std::size_t R, std::size_t C>
vector<T, R> column(matrix<T, R, C> const & m, std::size_t j)
{
vector<T, R> r;
for (std::size_t i = 0; i < R; ++i)
vector<T, R> r(m.coords.dimensions().rows);
for (std::size_t i = 0; i < m.rows(); ++i)
r[i] = m[i][j];
return r;
}
@ -373,8 +437,8 @@ namespace psemek::math
T frobenius_norm_sqr(matrix<T, R, C> const & m)
{
T r{0};
for (std::size_t i = 0; i < R; ++i)
for (std::size_t j = 0; j < C; ++j)
for (std::size_t i = 0; i < m.rows(); ++i)
for (std::size_t j = 0; j < m.columns(); ++j)
r += sqr(m[i][j]);
return r;
}
@ -389,12 +453,13 @@ namespace psemek::math
T trace(matrix<T, R, C> const & m)
{
T r{0};
for (std::size_t i = 0; i < std::min(R, C); ++i)
for (std::size_t i = 0; i < std::min(m.rows(), m.columns()); ++i)
r += m[i][i];
return r;
}
template <std::size_t R1, std::size_t C1, std::size_t R2, std::size_t C2, typename T, std::size_t R, std::size_t C>
requires (R1 != dynamic && C1 != dynamic && R2 != dynamic && C2 != dynamic)
matrix<T, R2 - R1, C2 - C1> submatrix(matrix<T, R, C> const & m)
{
static_assert(R1 < R2);
@ -403,13 +468,9 @@ namespace psemek::math
static_assert(C2 <= C);
matrix<T, R2 - R1, C2 - C1> result;
for (std::size_t i = 0; i < R2 - R1; ++i)
{
for (std::size_t j = 0; j < C2 - C1; ++j)
{
for (std::size_t i = 0; i < result.rows(); ++i)
for (std::size_t j = 0; j < result.columns(); ++j)
result[i][j] = m[i + R1][j + C1];
}
}
return result;
}
@ -426,12 +487,12 @@ namespace psemek::math
return r;
}
template <typename T, std::size_t N>
matrix<T, N, N> outer_product(vector<T, N> const & v1, vector<T, N> const & v2)
template <typename T, std::size_t R, std::size_t C>
matrix<T, R, C> outer_product(vector<T, R> const & v1, vector<T, C> const & v2)
{
matrix<T, N, N> r;
for (std::size_t i = 0; i < N; ++i)
for (std::size_t j = 0; j < N; ++j)
matrix<T, R, C> r({.rows = v1.coords.dimensions(), .columns = v2.coords.dimensions()});
for (std::size_t i = 0; i < r.rows(); ++i)
for (std::size_t j = 0; j < r.columns(); ++j)
r[i][j] = v1[i] * v2[j];
return r;
}
@ -439,9 +500,9 @@ namespace psemek::math
template <typename T, std::size_t R, std::size_t C>
std::ostream & operator << (std::ostream & os, matrix<T, R, C> const & m)
{
for (std::size_t i = 0; i < R; ++i)
for (std::size_t i = 0; i < m.rows(); ++i)
{
for (std::size_t j = 0; j < C; ++j)
for (std::size_t j = 0; j < m.columns(); ++j)
os << m[i][j] << ' ';
os << '\n';
}
@ -458,9 +519,9 @@ namespace psemek::math
template <typename T, std::size_t R, std::size_t C>
std::ostream & operator << (std::ostream & os, setw<T, R, C> const & w)
{
for (std::size_t i = 0; i < R; ++i)
for (std::size_t i = 0; i < w.m.rows(); ++i)
{
for (std::size_t j = 0; j < C; ++j)
for (std::size_t j = 0; j < w.m.columns(); ++j)
os << std::fixed << std::right << std::setw(w.w) << w.m[i][j] << ' ';
os << '\n';
}

View file

@ -19,17 +19,42 @@ namespace psemek::math
typename detail::array<T, N>::type coords;
point() = default;
using dynamic_dimensions_type = detail::dynamic_dimensions<N>;
// Internal
point(dynamic_dimensions_type dimensions)
: coords(dimensions)
{}
point()
: coords({})
{}
explicit point(std::size_t size) requires (N == dynamic)
: coords({.size = size})
{}
template <typename ... Args>
requires((sizeof...(Args) == N) && detail::all_convertible_to<T, Args...>::value)
requires(N != dynamic && sizeof...(Args) == N && detail::all_convertible_to<T, Args...>::value)
point(Args && ... args)
: coords{ static_cast<T>(std::forward<Args>(args))... }
{}
: coords({})
{
auto out = values().begin();
((*out++ = args), ...);
}
template <typename ... Args>
requires(N == dynamic && detail::all_convertible_to<T, Args...>::value)
point(Args && ... args)
: coords({.size = sizeof...(Args)})
{
auto out = values().begin();
((*out++ = args), ...);
}
std::size_t dimension() const
{
return N;
return coords.size;
}
T & operator[](std::size_t i)
@ -42,20 +67,56 @@ namespace psemek::math
return coords[i];
}
util::span<T> values()
{
return {&coords[0], dimension()};
}
util::span<T const> values() const
{
return {&coords[0], dimension()};
}
point copy() const
{
if constexpr (N == dynamic)
{
point result;
result.coords = coords.copy();
return result;
}
else
{
return *this;
}
}
point & operator += (vector<T, N> const & v);
point & operator -= (vector<T, N> const & v);
static point zero();
static point zero(dynamic_dimensions_type dimensions = {});
};
template <typename ... Args>
point(Args && ...) -> point<std::common_type_t<Args...>, sizeof...(Args)>;
template <typename T, std::size_t N>
point<T, N> point<T, N>::zero()
void check_dynamic_size(point<T, N> const & p1, point<T, N> const & p2)
{
point<T, N> p;
for (std::size_t i = 0; i < N; ++i)
check_dynamic_size(p1.coords.dimensions(), p2.coords.dimensions());
}
template <typename T, std::size_t N>
void check_dynamic_size(point<T, N> const & p, vector<T, N> const & v)
{
check_dynamic_size(p.coords.dimensions(), v.coords.dimensions());
}
template <typename T, std::size_t N>
point<T, N> point<T, N>::zero(dynamic_dimensions_type dimensions)
{
point<T, N> p(dimensions);
for (std::size_t i = 0; i < p.dimension(); ++i)
p[i] = T(0);
return p;
}
@ -63,8 +124,8 @@ namespace psemek::math
template <typename T1, typename T, std::size_t N>
point<T1, N> cast(point<T, N> const & p)
{
point<T1, N> r;
for (std::size_t i = 0; i < N; ++i)
point<T1, N> r(p.coords.dimensions());
for (std::size_t i = 0; i < r.dimension(); ++i)
r[i] = T1(p[i]);
return r;
}
@ -108,7 +169,9 @@ namespace psemek::math
template <typename T, std::size_t N>
std::strong_ordering operator <=> (point<T, N> const & p1, point<T, N> const & p2)
{
for (std::size_t i = 0; i < N; ++i)
check_dynamic_size(p1, p2);
for (std::size_t i = 0; i < p1.dimension(); ++i)
{
if (p1[i] < p2[i])
return std::strong_ordering::less;
@ -121,7 +184,8 @@ namespace psemek::math
template <typename T, std::size_t N>
point<T, N> operator + (point<T, N> const & p, vector<T, N> const & v)
{
point<T, N> r;
check_dynamic_size(p, v);
point<T, N> r(p.coords.dimensions());
for (std::size_t i = 0; i < N; ++i)
r[i] = p[i] + v[i];
return r;
@ -130,7 +194,8 @@ namespace psemek::math
template <typename T, std::size_t N>
point<T, N> operator + (vector<T, N> const & v, point<T, N> const & p)
{
point<T, N> r;
check_dynamic_size(p, v);
point<T, N> r(p.coords.dimensions());
for (std::size_t i = 0; i < N; ++i)
r[i] = v[i] + p[i];
return r;
@ -139,7 +204,8 @@ namespace psemek::math
template <typename T, std::size_t N>
point<T, N> operator - (point<T, N> const & p, vector<T, N> const & v)
{
point<T, N> r;
check_dynamic_size(p, v);
point<T, N> r(p.coords.dimensions());
for (std::size_t i = 0; i < N; ++i)
r[i] = p[i] - v[i];
return r;
@ -148,7 +214,8 @@ namespace psemek::math
template <typename T, std::size_t N>
vector<T, N> operator - (point<T, N> const & p1, point<T, N> const & p2)
{
vector<T, N> r;
check_dynamic_size(p1, p2);
vector<T, N> r(p1.coords.dimensions());
for (std::size_t i = 0; i < N; ++i)
r[i] = p1[i] - p2[i];
return r;
@ -199,8 +266,14 @@ namespace psemek::math
template <typename T, std::size_t N>
std::ostream & operator << (std::ostream & os, point<T, N> const & p)
{
if constexpr (N == 0)
{
os << "()";
return os;
}
os << '(' << p[0];
for (std::size_t i = 1; i < N; ++i)
for (std::size_t i = 1; i < p.dimension(); ++i)
os << ", " << p[i];
os << ')';
return os;
@ -209,7 +282,7 @@ namespace psemek::math
template <typename T, std::size_t N>
bool isfinite(point<T, N> const & p)
{
for (std::size_t i = 0; i < N; ++i)
for (std::size_t i = 0; i < p.dimension(); ++i)
if (!std::isfinite(p[i]))
return false;
return true;

View file

@ -2,11 +2,10 @@
#include <psemek/math/detail/array.hpp>
#include <psemek/util/hash.hpp>
#include <psemek/util/assert.hpp>
#include <psemek/util/span.hpp>
#include <iostream>
#include <cstddef>
#include <utility>
#include <type_traits>
#include <cmath>
#include <compare>
@ -22,33 +21,62 @@ namespace psemek::math
typename detail::array<T, N>::type coords;
vector() = default;
using dynamic_dimensions_type = detail::dynamic_dimensions<N>;
explicit vector(std::size_t size) requires (N != dynamic)
{
assert(size == N);
}
// Internal
vector(dynamic_dimensions_type dimensions)
: coords(dimensions)
{}
vector()
: coords({})
{}
explicit vector(std::size_t size) requires (N == dynamic)
: coords(size, T{})
: coords({.size = size})
{}
template <typename ... Args>
requires((sizeof...(Args) == N || N == dynamic) && detail::all_convertible_to<T, Args...>::value)
requires(N != dynamic && sizeof...(Args) == N && detail::all_convertible_to<T, Args...>::value)
vector(Args && ... args)
: coords{ static_cast<T>(std::forward<Args>(args))... }
{}
: coords({})
{
auto out = values().begin();
((*out++ = args), ...);
}
template <typename ... Args>
requires(N == dynamic && detail::all_convertible_to<T, Args...>::value)
vector(Args && ... args)
: coords({.size = sizeof...(Args)})
{
auto out = values().begin();
((*out++ = args), ...);
}
std::size_t dimension() const
{
if constexpr (N == dynamic)
{
return coords.size();
return coords.size;
}
else
T & operator[](std::size_t i)
{
return N;
return coords[i];
}
T const & operator[](std::size_t i) const
{
return coords[i];
}
util::span<T> values()
{
return {&coords[0], dimension()};
}
util::span<T const> values() const
{
return {&coords[0], dimension()};
}
vector copy() const
@ -65,33 +93,28 @@ namespace psemek::math
}
}
T & operator[](std::size_t i)
{
return coords[i];
}
T const & operator[](std::size_t i) const
{
return coords[i];
}
vector & operator *= (T const & s);
vector & operator /= (T const & s);
vector & operator += (vector const & v);
vector & operator -= (vector const & v);
static vector zero() requires (N != dynamic);
static vector zero(std::size_t size);
static vector zero(dynamic_dimensions_type dimensions = {});
};
template <typename ... Args>
vector(Args && ...) -> vector<std::common_type_t<Args...>, sizeof...(Args)>;
template <typename T, std::size_t N>
void check_dynamic_size(vector<T, N> const & v1, vector<T, N> const & v2)
{
check_dynamic_size(v1.coords.dimensions(), v2.coords.dimensions());
}
template <typename T1, typename T, std::size_t N>
vector<T1, N> cast(vector<T, N> const & v)
{
vector<T1, N> r(v.dimension());
vector<T1, N> r(v.coords.dimensions());
for (std::size_t i = 0; i < v.dimension(); ++i)
r[i] = static_cast<T1>(v[i]);
return r;
@ -100,7 +123,7 @@ namespace psemek::math
template <typename T, std::size_t N>
std::strong_ordering operator <=> (vector<T, N> const & v1, vector<T, N> const & v2)
{
check_dynamic_size(v1.dimension(), v2.dimension());
check_dynamic_size(v1, v2);
for (std::size_t i = 0; i < v1.dimension(); ++i)
{
@ -151,7 +174,7 @@ namespace psemek::math
template <typename T, std::size_t N>
vector<T, N> operator * (vector<T, N> const & v, T const & s)
{
vector<T, N> r(v.dimension());
vector<T, N> r(v.coords.dimensions());
for (std::size_t i = 0; i < v.dimension(); ++i)
r[i] = v[i] * s;
return r;
@ -160,7 +183,7 @@ namespace psemek::math
template <typename T, std::size_t N>
vector<T, N> operator * (T const & s, vector<T, N> const & v)
{
vector<T, N> r(v.dimension());
vector<T, N> r(v.coords.dimensions());
for (std::size_t i = 0; i < v.dimension(); ++i)
r[i] = s * v[i];
return r;
@ -169,7 +192,7 @@ namespace psemek::math
template <typename T, std::size_t N>
vector<T, N> operator / (vector<T, N> const & v, T const & s)
{
vector<T, N> r(v.dimension());
vector<T, N> r(v.coords.dimensions());
for (std::size_t i = 0; i < v.dimension(); ++i)
r[i] = v[i] / s;
return r;
@ -178,7 +201,7 @@ namespace psemek::math
template <typename T, std::size_t N>
vector<T, N> operator - (vector<T, N> const & v)
{
vector<T, N> r(v.dimension());
vector<T, N> r(v.coords.dimensions());
for (std::size_t i = 0; i < v.dimension(); ++i)
r[i] = -v[i];
return r;
@ -187,9 +210,9 @@ namespace psemek::math
template <typename T, std::size_t N>
vector<T, N> operator + (vector<T, N> const & v1, vector<T, N> const & v2)
{
check_dynamic_size(v1.dimension(), v2.dimension());
check_dynamic_size(v1, v2);
vector<T, N> r(v1.dimension());
vector<T, N> r(v1.coords.dimensions());
for (std::size_t i = 0; i < v1.dimension(); ++i)
r[i] = v1[i] + v2[i];
return r;
@ -198,9 +221,9 @@ namespace psemek::math
template <typename T, std::size_t N>
vector<T, N> operator - (vector<T, N> const & v1, vector<T, N> const & v2)
{
check_dynamic_size(v1.dimension(), v2.dimension());
check_dynamic_size(v1, v2);
vector<T, N> r(v1.dimension());
vector<T, N> r(v1.coords.dimensions());
for (std::size_t i = 0; i < v1.dimension(); ++i)
r[i] = v1[i] - v2[i];
return r;
@ -231,24 +254,10 @@ namespace psemek::math
}
template <typename T, std::size_t N>
vector<T, N> vector<T, N>::zero() requires (N != dynamic)
vector<T, N> vector<T, N>::zero(dynamic_dimensions_type dimensions)
{
vector<T, N> r;
for (std::size_t i = 0; i < N; ++i)
r[i] = T(0);
return r;
}
template <typename T, std::size_t N>
vector<T, N> vector<T, N>::zero(std::size_t size)
{
if constexpr (N != dynamic)
{
assert(N == size);
}
vector<T, N> r(size);
for (std::size_t i = 0; i < size; ++i)
vector<T, N> r(dimensions);
for (std::size_t i = 0; i < r.dimension(); ++i)
r[i] = T(0);
return r;
}
@ -256,7 +265,7 @@ namespace psemek::math
template <typename T, std::size_t N>
T dot(vector<T, N> const & v1, vector<T, N> const & v2)
{
check_dynamic_size(v1.dimension(), v2.dimension());
check_dynamic_size(v1, v2);
T r{};
for (std::size_t i = 0; i < v1.dimension(); ++i)
@ -332,8 +341,8 @@ namespace psemek::math
template <typename T, std::size_t N>
vector<T, N> ort(vector<T, N> const & v)
{
vector<T, N> result = vector<T, N>::zero(v.dimension());
if constexpr ((N % 2) == 0)
vector<T, N> result = vector<T, N>::zero(v.coords.dimensions());
if ((N != dynamic && (N % 2) == 0) || (N == dynamic && (v.dimension() % 2) == 0))
{
for (std::size_t i = 0; i < v.dimension(); i += 2)
{
@ -414,9 +423,9 @@ namespace psemek::math
{
using R = decltype(f(v0[0], vs[0]...));
(check_dynamic_size(v0.dimension(), vs.dimension()), ...);
(check_dynamic_size(v0, vs), ...);
vector<R, N> result(v0.dimension());
vector<R, N> result(v0.coords.dimensions());
for (std::size_t i = 0; i < v0.dimension(); ++i)
result[i] = f(v0[i], vs[i]...);
return result;