chsvmath\algebraic_traits.cpp

Demonstrates usage of algebraic traits templating to select a suitable form of an algebraic algorithm to provide the best performance.

#include "chsvmath.h"
#include <cstdint>
#include <iostream>
#include <utility>
#include <type_traits>
 
template <class T>
constexpr auto AdditiveInverse(const T& val) noexcept
{
    return -val;
}
 
template <class T>
auto MultiplicativeInverse(const T& val) noexcept -> decltype(std::declval<const T&>().MultiplicativeInverse())
{
    return val.MultiplicativeInverse();
}
 
template <class T>
auto MultiplicativeInverse(const T& val) noexcept -> std::enable_if_t<std::is_arithmetic_v<T>, T>
{
    return T(1) / val;
}
 
template <class T, std::size_t nRows, std::size_t nCols>
T Determinant(const T (&Matrix)[nRows][nCols], Chusov::Math::rng_algebraic_tag)
{
    //For Rng's a division operation is not defined, hence Gaussian transformation of matrix rows is not applicable,
    //and only complex O(n!) algorithm can be used to calculate the determinant.
    static_assert(nRows == nCols, "The matrix must be square.");
    T result = T();
    T Minor[nRows - 1][nCols - 1];
    for (std::size_t i = 0; i < nCols; ++i)
    {
        for (std::size_t iRow = 1; iRow < nRows; ++iRow)
            std::copy_if(&Matrix[iRow][0], &Matrix[iRow][nCols], &Minor[iRow - 1][0], [iRow, i, &Matrix](const T& val) {return &val - &Matrix[iRow][0] != i;});
        result += ((i & 1) == 0?Determinant(std::move(Minor)):AdditiveInverse(Determinant(std::move(Minor)))) * Matrix[0][i];
    }
    return result;
}
 
template <class T>
T Determinant(const T (&Matrix)[1][1], Chusov::Math::rng_algebraic_tag)
{
    return Matrix[0][0];
}
 
template <class T, std::size_t nRows, std::size_t nCols>
T Determinant(T (&&Matrix)[nRows][nCols], Chusov::Math::division_ring_algebraic_tag)
{
    static_assert(nRows == nCols, "The matrix has to be square.");
    //Division rings define inversion over multiplication which gives the opportunity to apply Gaussian transformation to the matrix and hence reduce
    //the complexity of the algorithm to O(n^3)
    T result = 1;
    bool fNeg = false;
    for (std::size_t iCol = 0; iCol < nCols; ++iCol)
    {
        auto iRow = iCol;
        if (Matrix[iRow][iCol] == T())
        {
            auto nzi = std::find_if(&Matrix[iRow + 1], &Matrix[nRows], [iCol](const T (&val)[nCols]) {return val[iCol] != T();}) - &Matrix[0];
            if (nzi == nRows)
                return T();
            std::swap(Matrix[iRow], Matrix[nzi]);
            iRow = nzi;
            fNeg = !fNeg;
        }
        for (++iRow; iRow < nRows; ++iRow)
        {
            /*
            For the following rows:
            [a b ... c]
            [d e ... f]
            The following transformation is performed during the iteration:
            [a                b              ... c             ]
            [d-(d*(a^-1))*a   e-(d*(a^-1))*b ... f-(d*(a^-1))*c]
            */
            T mult = Matrix[iRow][iCol] * MultiplicativeInverse(Matrix[iCol][iCol]);
            for (auto iRowCol = iCol; iRowCol < nRows; ++iRowCol)
                Matrix[iRow][iRowCol] = Matrix[iRow][iRowCol] - mult * Matrix[iCol][iRowCol];
        }
        result *= Matrix[iCol][iCol];
    }
    return fNeg?AdditiveInverse(result):result;
}
 
template <class T, std::size_t nRows, std::size_t nCols>
T Determinant(const T (&Matrix)[nRows][nCols], Chusov::Math::division_ring_algebraic_tag tag)
{
    T ModifiableCopy[nRows][nCols];
    std::copy(&Matrix[0][0], &Matrix[nRows][0], &ModifiableCopy[0][0]);
    return Determinant(std::move(ModifiableCopy), tag);
}
 
template <class T, std::size_t nRows, std::size_t nCols>
T Determinant(const T (&Matrix)[nRows][nCols])
{
    return Determinant(Matrix, typename Chusov::Math::algebraic_traits<T>::algebraic_category());
}
 
class prime_field_11;
 
template <class _Char_t, class _Traits_t>
std::basic_ostream<_Char_t, _Traits_t>& operator<<(std::basic_ostream<_Char_t, _Traits_t>& os, const prime_field_11& x);
 
class prime_field_11 //finite field
{
    static const std::size_t _size = 11;
    static const unsigned values[_size];
    static const unsigned additive_inversion[_size];
    static const unsigned multiplicative_inversion[_size];
    unsigned m_val;
    explicit operator unsigned() const {return m_val;}
public:
    typedef Chusov::Math::field_algebraic_tag algebraic_category; //algebraic traits tag
    prime_field_11()=default;
    /*
    During a prime_field_11(Chusov::Math::additive_identity) call the Chusov::Math::additive_identity value implicitly converts to 0.
    During a prime_field_11(Chusov::Math::multiplicative_identity) call the Chusov::Math::multiplicative_identity value implicitly converts to 1.
    */
    prime_field_11(unsigned val):m_val(val % _size){}
    
    //mandatory methods for the class to adhere the field_concept:
    prime_field_11& operator+=(const prime_field_11& addend) {(m_val += addend.m_val) %= _size; return *this;}
    prime_field_11 operator+(const prime_field_11& addend) const {return (m_val + addend.m_val) % _size;}
    prime_field_11& operator-=(const prime_field_11& subtrahend)
    {
        m_val = subtrahend.m_val <= m_val?m_val - subtrahend.m_val:additive_inversion[subtrahend.m_val - m_val];
        return *this;
    }
    prime_field_11 operator-(const prime_field_11& subtrahend) const
    {
        return subtrahend.m_val <= m_val?m_val - subtrahend.m_val:additive_inversion[subtrahend.m_val - m_val];
    }
    prime_field_11 operator-() const {return additive_inversion[m_val];}
    prime_field_11& operator*=(const prime_field_11& multiplier) {(m_val *= multiplier.m_val) %= _size; return *this;}
    prime_field_11 operator*(const prime_field_11& multiplier) const {return (m_val * multiplier.m_val) % _size;}
    prime_field_11& operator/=(const prime_field_11& dividend) {return *this *= dividend.MultiplicativeInverse();}
    prime_field_11 operator/(const prime_field_11& dividend) const {return *this * dividend.MultiplicativeInverse();}
    prime_field_11 MultiplicativeInverse() const {return multiplicative_inversion[m_val];}
    prime_field_11& MultiplicativeInverseMe() {m_val = multiplicative_inversion[m_val]; return *this;}
    bool operator==(const prime_field_11& _right) const {return m_val == _right.m_val;}
    template <class _Char_t, class _Traits_t>
    friend std::basic_ostream<_Char_t, _Traits_t>& operator<<(std::basic_ostream<_Char_t, _Traits_t>& os, const prime_field_11& x);
};
 
template <class _Char_t, class _Traits_t>
std::basic_ostream<_Char_t, _Traits_t>& operator<<(std::basic_ostream<_Char_t, _Traits_t>& os, const prime_field_11& x)
{
    os << unsigned(x);
    return os;
}
 
const unsigned prime_field_11::values[prime_field_11::_size]                    = {            0,  1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
const unsigned prime_field_11::additive_inversion[prime_field_11::_size]        = {            0, 10, 9, 8, 7, 6, 5, 4, 3, 2,  1};
const unsigned prime_field_11::multiplicative_inversion[prime_field_11::_size]  = {(unsigned) -1,  1, 6, 4, 3, 9, 2, 8, 7, 5, 10};
 
 
class integral_ring_12;
 
template <class _Char_t, class _Traits_t>
std::basic_ostream<_Char_t, _Traits_t>& operator<<(std::basic_ostream<_Char_t, _Traits_t>& os, const integral_ring_12& x);
 
class integral_ring_12
{
    static constexpr std::size_t _size = 12;
    unsigned m_val;
    explicit operator unsigned() const {return m_val;}
public:
    typedef Chusov::Math::commutative_ring_algebraic_tag algebraic_category; //algebraic traits tag
    integral_ring_12()=default;
    /*
    During a integral_ring_12(Chusov::Math::additive_identity) call the Chusov::Math::additive_identity value implicitly converts to 0.
    During a integral_ring_12(Chusov::Math::multiplicative_identity) call the Chusov::Math::multiplicative_identity value implicitly converts to 1.
    */
    integral_ring_12(unsigned val):m_val(val % _size){}
    
    //mandatory methods for the class to adhere the ring_concept:
    integral_ring_12& operator+=(const integral_ring_12& addend) {(m_val += addend.m_val) %= _size; return *this;}
    integral_ring_12 operator+(const integral_ring_12& addend) const {return (m_val + addend.m_val) % _size;}
    integral_ring_12& operator-=(const integral_ring_12& subtrahend)
    {
        m_val = subtrahend.m_val <= m_val?m_val - subtrahend.m_val:(unsigned) _size - (subtrahend.m_val - m_val);
        return *this;
    }
    integral_ring_12 operator-(const integral_ring_12& subtrahend) const
    {
        return subtrahend.m_val <= m_val?m_val - subtrahend.m_val:(unsigned) _size - (subtrahend.m_val - m_val);
    }
    integral_ring_12 operator-() const {return m_val == 0?0:(unsigned) _size - m_val;}
    integral_ring_12& operator*=(const integral_ring_12& multiplier) {(m_val *= multiplier.m_val) %= _size; return *this;}
    integral_ring_12 operator*(const integral_ring_12& multiplier) const {return (m_val * multiplier.m_val) % _size;}
    bool operator==(const integral_ring_12& _right) const {return m_val == _right.m_val;}
    template <class _Char_t, class _Traits_t>
    friend std::basic_ostream<_Char_t, _Traits_t>& operator<<(std::basic_ostream<_Char_t, _Traits_t>& os, const integral_ring_12& x);
};
 
template <class _Char_t, class _Traits_t>
std::basic_ostream<_Char_t, _Traits_t>& operator<<(std::basic_ostream<_Char_t, _Traits_t>& os, const integral_ring_12& x)
{
    os << unsigned(x);
    return os;
}
 
namespace std
{
    template <class container_t>
    std::string to_string(const container_t& cont)
    {
        std::ostringstream os;
        os << '{';
        for (const auto& val:cont)
        {
            if (val == *std::begin(cont))
                os << val;
            else
                os << ", " << val;
        }
        os << '}';
        return os.str();
    }
}
 
int main(int argc, char** argv)
{
    //Calls for the first overload of the Determinant function
    integral_ring_12 rnz[3][3] = 
    {
        {1, 2, 3},
        {4, 3, 2},
        {7, 1, 5}
    };
    integral_ring_12 rz[3][3] =
    {
        {1, 2, 3},
        {4, 5, 6},
        {7, 8, 9}
    };
    std::cout << "Determinant of\n";
    for (auto& line:rnz)
        std::cout << std::to_string(line) << std::endl;
    std::cout << "is " << Determinant(rnz) << std::endl; //Calls the first overload of the Determinant function expected result is 10
    
    std::cout << "Determinant of\n";
    for (auto& line:rz)
        std::cout << std::to_string(line) << std::endl;
    std::cout << "is " << Determinant(rz) << std::endl;  //Calls the first overload of the Determinant function expected result is 0
    
    //Calls for the second overload of the Determinant function
    prime_field_11 fnz[3][3] = 
    {
        {1, 2, 3},
        {4, 3, 2},
        {7, 1, 5}
    };
    prime_field_11 fz[3][3] =
    {
        {1, 2, 3},
        {4, 5, 6},
        {7, 8, 9}
    };
    std::cout << "Determinant of\n";
    for (auto& line:fnz)
        std::cout << std::to_string(line) << std::endl;
    std::cout << "is " << Determinant(fnz) << std::endl; //Calls the second overload of the Determinant function expected result is 5
    
    std::cout << "Determinant of\n";
    for (auto& line:fz)
        std::cout << std::to_string(line) << std::endl;
    std::cout << "is " << Determinant(fz) << std::endl;  //Calls the second overload of the Determinant function expected result is 0
    
    //Calls for the second overload of the Determinant function
    double dnz[3][3] = 
    {
        {1, 2, 3},
        {4, 3, 2},
        {7, 1, 5}
    };
    double dz[3][3] =
    {
        {1, 2, 3},
        {4, 5 ,6},
        {7, 8, 9}
    };
    std::cout << "Determinant of\n";
    for (auto& line:dnz)
        std::cout << std::to_string(line) << std::endl;
    std::cout << "is " << Determinant(dnz) << std::endl; //Calls the first overload of the Determinant function expected result is -50
    
    std::cout << "Determinant of\n";
    for (auto& line:dz)
        std::cout << std::to_string(line) << std::endl;
    std::cout << "is " << Determinant(dz) << std::endl;  //Calls the first overload of the Determinant function expected result is 0
    
    return 0;
}