algebra/include/nil/crypto3/algebra/matrix/math.hpp

Go to the documentation of this file.

//---------------------------------------------------------------------------//
// Copyright (c) 2020-2021 Mikhail Komarov <nemo@nil.foundation>
// Copyright (c) 2020-2021 Nikita Kaskov <nbering@nil.foundation>
// Copyright (c) 2020-2021 Ilias Khairullin <ilias@nil.foundation>
//
// MIT License
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
//---------------------------------------------------------------------------//
 
#ifndef CRYPTO3_ALGEBRA_MATRIX_MATH_HPP
#define CRYPTO3_ALGEBRA_MATRIX_MATH_HPP
 
#include <algorithm>
 
#include <nil/crypto3/detail/assert.hpp>
 
#include <nil/crypto3/algebra/vector/vector.hpp>
#include <nil/crypto3/algebra/vector/math.hpp>
 
#include <nil/crypto3/algebra/matrix/matrix.hpp>
#include <nil/crypto3/algebra/matrix/utility.hpp>
 
#include <nil/crypto3/algebra/scalar/math.hpp>
 
namespace nil {
    namespace crypto3 {
        namespace algebra {
            template<typename T, std::size_t M, std::size_t N>
            constexpr matrix<T, M, N> conj(const matrix<T, M, N> &m) {
                return elementwise(algebra::conj<T>, m);
            }
 
            template<typename T, std::size_t M, std::size_t N>
            constexpr matrix<nil::crypto3::detail::remove_complex_t<T>, M, N> real(const matrix<T, M, N> &m) {
                return elementwise([](auto i) { return std::real(i); }, m);
            }
 
            template<typename T, std::size_t M, std::size_t N>
            constexpr matrix<nil::crypto3::detail::remove_complex_t<T>, M, N> imag(const matrix<T, M, N> &m) {
                return elementwise([](auto i) { return std::imag(i); }, m);
            }
 
            template<typename T, std::size_t M, std::size_t N>
            constexpr matrix<T, N, M> transpose(const matrix<T, M, N> &m) {
                return generate<N, M>([&m](auto i, auto j) { return m[j][i]; });
            }
 
            template<typename T, std::size_t M, std::size_t N>
            constexpr matrix<T, N, M> hermitian(const matrix<T, M, N> &m) {
                return transpose(conj(m));
            }
 
            template<typename T, std::size_t M, std::size_t N, std::size_t P>
            constexpr matrix<T, M, P> matmul(const matrix<T, M, N> &a, const matrix<T, N, P> &b) {
                return generate<M, P>([&a, &b](auto i, auto j) { return algebra::sum(a.row(i) * b.column(j)); });
            }
 
            template<typename T, std::size_t M, std::size_t N>
            constexpr vector<T, N> vectmatmul(const vector<T, M> &v, const matrix<T, M, N> &m) {
                return generate<N>([&v, &m](auto i) { return sum(v * m.column(i)); });
            }
 
            template<typename T, std::size_t M, std::size_t N>
            constexpr vector<T, M> matvectmul(const matrix<T, M, N> &m, const vector<T, N> &v) {
                return generate<M>([&v, &m](auto i) { return sum(m.row(i) * v); });
            }
 
            template<typename T, std::size_t M, std::size_t N, std::size_t P, std::size_t Q>
            constexpr matrix<T, M * P, N * Q> kron(const matrix<T, M, N> &a, const matrix<T, P, Q> &b) {
                return generate<M * P, N * Q>([&a, &b](auto i, auto j) { return a[i / P][j / Q] * b[i % P][j % Q]; });
            }
 
            template<typename T, std::size_t M, std::size_t N>
            constexpr T macs(const matrix<T, M, N> &m) {
                return max(generate<N>([&m](std::size_t i) { return sum(abs(m.column(i))); }));
            }
 
            template<typename T, std::size_t M, std::size_t N>
            constexpr T mars(const matrix<T, M, N> &m) {
                return max(generate<M>([&m](std::size_t i) { return sum(abs(m.row(i))); }));
            }
 
            template<typename T, std::size_t M, std::size_t N>
            constexpr std::tuple<matrix<T, M, N>, std::size_t, T> gauss_jordan_impl(matrix<T, M, N> m, T tolerance) {
                // CRYPTO3_DETAIL_ASSERT_FLOATING_POINT(T)
                // CRYPTO3_DETAIL_ASSERT_REAL(T)
 
                // Define function for determining if an element is negligible
                auto negligible = [&tolerance](const T &v) { return abs(v) < tolerance; };
 
                T det = 1;
                std::size_t rank = 0;
                std::size_t i = 0, j = 0;
                while (i < M && j < N) {
                    // Choose largest magnitude as pivot to avoid adding different magnitudes
                    for (std::size_t ip = i + 1; ip < M; ++ip) {
                        if (abs(m[ip][j]) > abs(m[i][j])) {
                            for (std::size_t jp = 0; jp < N; ++jp) {
                                auto tmp = m[ip][jp];
                                m[ip][jp] = m[i][jp];
                                m[i][jp] = tmp;
                            }
                            det = -det;
                            break;
                        }
                    }
 
                    // If m_ij is still 0, continue to the next column
                    if (!negligible(m[i][j])) {
                        // Scale m_ij to 1
                        auto s = m[i][j];
                        for (std::size_t jp = 0; jp < N; ++jp)
                            m[i][jp] /= s;
                        det /= s;
 
                        // Eliminate other values in the column
                        for (std::size_t ip = 0; ip < M; ++ip) {
                            if (ip == i)
                                continue;
                            if (!negligible(m[ip][j])) {
                                auto s = m[ip][j];
                                [&]() {    // wrap this in a lambda to get around a gcc bug
                                    for (std::size_t jp = 0; jp < N; ++jp)
                                        m[ip][jp] -= s * m[i][jp];
                                }();
                            }
                        }
 
                        // Increment rank
                        ++rank;
 
                        // Select next row
                        ++i;
                    }
                    ++j;
                }
                det = (rank == M) ? det : 0;
                return {m, rank, det};
            }
 
            template<typename T, std::size_t M, std::size_t N>
            constexpr std::tuple<matrix<T, M, N>, std::size_t, T> gauss_jordan_impl(const matrix<T, M, N> &m) {
                T tol = T(std::max(N, M)) * std::numeric_limits<T>::epsilon() * mars(m);
                return gauss_jordan_impl(m, tol);
            }
 
            template<typename T, std::size_t M, std::size_t N>
            constexpr matrix<T, M, N> rref(const matrix<T, M, N> &m) {
                return std::get<0>(gauss_jordan_impl(m));
            }
 
            template<typename T, std::size_t M, std::size_t N>
            constexpr matrix<T, M, N> rref(const matrix<T, M, N> &m, T tolerance) {
                return std::get<0>(gauss_jordan_impl(m), tolerance);
            }
 
            template<typename T, std::size_t M, std::size_t N>
            constexpr std::size_t rank(const matrix<T, M, N> &m) {
                return std::get<1>(gauss_jordan_impl(m));
            }
 
            template<typename T, std::size_t M>
            constexpr T det(const matrix<T, M, M> &m) {
                return std::get<2>(gauss_jordan_impl(m));
            }
 
            template<typename T, std::size_t M>
            constexpr matrix<T, M, M> inverse(const matrix<T, M, M> &m) {
                if (rank(m) < M)
                    throw "matrix is not invertible";
                return submat<M, M>(rref(horzcat(m, get_identity<T, M>())), 0, M);
            }
 
            template<typename T, std::size_t M>
            constexpr T trace(const matrix<T, M, M> &m) {
                return sum(generate<M>([&m](std::size_t i) { return m[i][i]; }));
            }
 
        }    // namespace algebra
    }        // namespace crypto3
}    // namespace nil
 
#endif    // CRYPTO3_ALGEBRA_MATRIX_MATH_HPP