Armadillo C++ : Linear Combination with modulus ca

I want to extract linear combinations from matrices but by performing combinations in modulus.

Let us consider the calculation modulus 5, we then have the following for the addition:

 + |  0 1 2 3 4
 --+-----------
 0 | 0 1 2 3 4
 1 | 1 2 3 4 0
 2 | 2 3 4 0 1
 3 | 3 4 0 1 2 
 4 | 4 0 1 2 3

and this table for the multiplication:

 * | 0 1 2 3 4
 --+-----------
 0 | 0 0 0 0 0
 1 | 0 1 2 3 4
 2 | 0 2 4 1 3
 3 | 0 3 1 4 2
 4 | 0 4 3 2 1

So let us take an example: Let us consider the following matrix:

E = 2 1 3 2 0
    4 3 0 1 1

Then we can obtain the triangulation matrix by applying a LU decomposition (https://en.wikipedia.org/wiki/LU_decomposition) (or a Gaussian Elimination), which is the following:

T = 1 0 0 0 0 
    2 1 0 0 0

and finally, the matrix that I want to extract, which is the one storing the linear combinations:

CL = 3 2 3 0 3
     0 1 1 3 4
     0 0 1 0 0 
     0 0 0 1 0
     0 0 0 0 1

So basically, the algorithm should work as follows:

 Input: a matrix E with n rows and m columns, and p, a prime number.

 * We perform a Gaussian elimination/LU decomposition to obtain the lower-triangulation matrix T. 
   But all the calculus are done modulo 'p'. 

 Output: T (with the same size as E, n rows m columns). 
         CL (with a size m rows, m columns), 
         which is basically the identity matrix on which we 
         applied all the modifications that were performed on E to obtain T.

Alright, so now we have the context, let me explain the problem. I started to do it using the Armadillo library (http://arma.sourceforge.net/), but I did not find any solution on the library to perform the calculus on a mathematical Field p. I easily found the LU method to obtain the lower-triangle matrix, but the calculations are performed in the real.

#include <iostream>
#include <armadillo>

using namespace arma;
using namespace std;

int main(int argc,char** argv)
{

    mat A = mat({{2,1,3,2,0},{4,3,0,1,1}});

    mat L, U, P;

    lu(L, U, P, A);

    cout << L << endl;

    return 0;
}

With the following, you obtain the lower-triangle matrix 'L' but in the real calculus. Thus you obtain:

 T' =  1   0
       1/2 1

Is there any technique to perform the computation in a modulus way?

EDIT The Armadillo library is not able to do it. I developed my own LU decomposition in modulus but there is still a bug there. I asked a new question here Linear Combination C++ in modulus, hoping to solve it.

First of all: drop the using namespaces, code can become completely unreadable if you do that.

I haven't used Armadillo yet. But I have looked at the documentation, and it seems templated to me.

Now things are getting a bit wild. The type you use, arma::mat seems to be a typedef on arma::Mat<double>.

The high-level function arma::lu isn't properly documented. It obviously does an LU-decomposition, but I don't know if the function is templated. If it is, i.e., you cannot just call it with double mats but also other types, you might have a shot using a custom type representing the field (since 5 is prime, otherwise you'd be completely lost) of calculations modulo 5. Meaning you write a class, let's call it IntMod5 and define all required operators for this class, meaning all operators that IntMod5 uses. For example, you'd need to define operator/(), e.g. by making a table of inverses of 4 of the 5 elements of the field (0 has none), i.e. 1->1, 2->3, 3->2, 4->4, and define

IntMod5 operator/(const IntMod5& o) const
{
    return IntMod5((this->value*inverse(o.value))%5);
}

This is just one example, you likely need to define all arithmetic operators, binary and unary, and possibly more such as comparison (LU decomposition might use finding good pivot elements). If you're lucky and the library is written in a way that it works for any field, not just floating point, you have a chance.

Before you go through all the work, you should use a trivial class simply wrapping double and check if arma::Mat or arma::lu do any type checks blocking you out.

If either of these fails, you'll likely have to write your own LU decomposition modulo 5 or find another library that supports it.