To work with integers (of arbitrary precision), use the following:

#include <lela/ring/integers.h>

using namespace LELA;

int main ()
{
	Integers Z;
	...
	return 0;
}

The element-type of Integers is called integer. It is defined in lela/integer.h.

To work with rational numbers, use the following:

#include <lela/ring/rationals.h>

using namespace LELA;

int main ()
{
	Rationals F;
	...
	return 0;
}

The element-type of Rationals is called RationalElement. It is defined in lela/element/rational.h.

Here is an example of the declaration of the field of integers modulo 101:

#include <lela/integer.h>
#include <lela/ring/modular.h>

using namespace LELA;

int main ()
{
	Modular<uint8> F (101);
	...
	return 0;
}

The type Modular is parametrised by the element-type, which may in principle be any C++-type which supports basic arithmetic operations (i.e. binary +, -, *, /, and %). The parameter passed to the constructor is the modulus. It may not exceed the maximum value allowed by the type.

The choice of which element-type to use is based on the size of the modulus (since obviously some moduli are too big for some element-types) as well as whether one wishes to use an integral of a floating-point representation. A floating-point representation generally uses more memory for the same modulus but enables the use of highly tuned BLAS-implementations which can greatly improve performance as compared with the integral representation.

While with a given element-type one can use any modulus up to the maximum (see table below), we recommend using a somewhat larger element-type than is necessary, since this can improve performance greatly. Here we present some recommendations.

The definitions of the types above (other than float and double) are in lela/integer.h. The header-file lela/ring/modular.h declares the class Modular.

The following code permits working over GF2:

#include <lela/ring/gf2.h>

using namespace LELA;

int main ()
{
	GF2 F;
	...
	return 0;
}

It is also possible to work over GF(2) using the Modular class described above, but the specialised class GF2 is vastly more efficient. Its use is highly recommended.

The element-type of GF2 is bool. The value true corresponds to the field-element 1 and false to 0.

If desired, it is possible to wrap a C++-type (builtin or user-defined) which supports standard arithmetic into LELA-ring so that it may be used within LELA. For example, wrapping float or double permits that LELA be used for numerical linear algebra. To do this, use the following;

#include <lela/ring/unparametric.h>

using namespace LELA;

int main ()
{
	TypeWrapperRing<float> R;
	...
	return 0;
}

The class TypeWrapperRing is parametrised by the type being wrapped, which is then the element-type of the resulting ring.

Once a ring is declared, it may be used to perform arithmetic:

#include <iostream>
#include <lela/integer.h>
#include <lela/ring/modular.h>

using namespace LELA;

int main ()
{
	Modular<uint8> F (101);

	uint8 a, b, c, d; // Declare ring-elements
	F.init (a, -1);   // Set a to -1 in the ring
	F.init (b, 1);    // Set b to 1 in the ring
	F.add (c, a, b);  // Set c <- a + b
	F.sub (c, a, b);  // Set c <- a - b
	F.mul (c, a, b);  // Set c <- a * b
	F.neg (c, b);     // Set c <- -b

	F.axpy (d, a, b, c);  // Set d <- a * b + c

	if (F.inv (c, a)) // Is a invertible in the ring? If so, set c to the inverse
		std::cout << "a is invertible" << std::endl;
	else
		std::cout << "a is not invertible" << std::endl;

	if (F.div (c, a, b)) // Does the quotient a / b exist in the ring? If so, set c to the quotient
		std::cout << "a / b exists" << std::endl;
	else
		std::cout << "a / b does not exist" << std::endl;

	return 0;
}

All arithmetic operations take a reference to the element in which to write the result as their first parameter. As compared with the seemingly more natural approach of returning the result, this makes things more efficient in the case that the element-type is complicated (such as a GMP-integer).

The method init sets its first argument to the ring-element which is the image of its second argument under the canonical map from the integers to the ring. The second argument may be of any integral type.

All arithmetic-operations have in-place variants:

#include <lela/integer.h>
#include <lela/ring/modular.h>

using namespace LELA;

int main ()
{
	Modular<uint8> F (101);

	uint8 a, b, c, d; // Declare ring-elements
	F.init (a, -1);   // Set a to -1 in the ring
	F.copy (b, a);    // Copy a to b
	F.negin (b);      // Replace b to -b
	F.init (c, 2);    // Set c to 2 in the ring
	F.addin (a, b);   // Set a <- a + b
	F.subin (a, b);   // Set a <- a - b
	F.mulin (a, b);   // Set a <- a * b

	F.axpyin (a, b, c);  // Set a <- a + b * c

	if (F.invin (a)) // Is a invertible in the ring? If so, replace it with its inverse
		std::cout << "a is invertible" << std::endl;
	else
		std::cout << "a is not invertible" << std::endl;

	if (F.divin (a, b)) // Does the quotient a / b exist in the ring? If so, replace a by a / b
		std::cout << "a / b exists" << std::endl;
	else
		std::cout << "a / b does not exist" << std::endl;

	return 0;
}

The method copy makes a copy of its second argument and stores it in its first argument. In generic code, it is recommended to use copy rather than the assignment-operator, since the latter may result in only a shallow copy, depending on the element-type. In the case of Modular<uint8> above, the call copy (b, a) is entirely equivalent to the assignment b = a.

Ring-elements may be read from or written to a C++-stream with the methods read and write:

#include <lela/integer.h>
#include <lela/ring/modular.h>

using namespace LELA;

int main ()
{
	Modular<uint8> F (101);

	uint8 a;

	std::cout << "Enter a value: ";
	F.read (std::cin, a);

	std::cout << "You entered ";
	F.write (std::cout, a) << std::endl;

	return 0;
}

The precise meaning of these methods depends on the ring. The method Modular::read interprets the input as an integer and maps it to the modular field under the canonical homomorphism.

Equality of elements may be checked with the method areEqual:

#include <lela/integer.h>
#include <lela/ring/modular.h>

using namespace LELA;

int main ()
{
	Modular<uint8> F (101);

	uint8 a, b;

	F.init (a, 50);

	std::cout << "Enter a value: ";
	F.read (std::cin, b);

	if (F.areEqual (a, b))
		std::cout << "You guessed correctly!" << std::endl;

	return 0;
}

Similarly, the methods isZero and isOne determine whether a ring-element is equal to zero or one, respectively, in the ring.

#include <lela/integer.h>
#include <lela/ring/modular.h>

using namespace LELA;

int main ()
{
	Modular<uint8> F (101);

	uint8 a;

	std::cout << "Enter a value: ";
	F.read (std::cin, a);

	if (F.isZero (a))
		std::cout << "You entered zero in the ring" << std::endl;
	if (F.isOne (a))
		std::cout << "You entered one in the ring" << std::endl;

	return 0;
}

The methods zero, one, and minusOne return read-only references to respectively the ring's 0, 1, and -1 elements.

#include <lela/integer.h>
#include <lela/ring/modular.h>

using namespace LELA;

int main ()
{
	Modular<uint8> F (101);

	uint8 a;

	std::cout << "Enter a value: ";
	F.read (std::cin, a);

	if (F.areEqual (a, F.one ()))
		std::cout << "You entered one in the ring" << std::endl;
	if (F.areEqual (a, F.zero ()))
		std::cout << "You entered zero in the ring" << std::endl;
	if (F.areEqual (a, F.minusOne ()))
		std::cout << "You entered minus one in the ring" << std::endl;

	return 0;
}

When writing algorithms, it is desireable that the code not depend on operating over a given ring. Thus the relevant function or class should be parametrised on the ring-type. Ring-elements may be declared with typename Ring::Element, as in the following:

template <class Ring>
void MyAlgorithm (const Ring &R, typename Ring::Element &v)
{
	typename Ring::Element a, b;

	R.init (a, 1);
	R.init (b, 103);
	R.add (v, a, b);
}

LELA supports five basic vector representation-types: two for general rings and three over GF(2). Any C++-type which meets the required interface of one of these types can be used in any container or algorithm which accepts that type.

#include <vector>

template <class Ring>
void MyAlgorithm (const Ring &R, size_t n)
{
	std::vector<typename Ring::Element> v (n); // A dense vector in R^n

	// Initialise some entries
	R.init (v[0], 1);
	R.init (v[1], -1);
	R.init (v[2], -3);
}

#include <lela/vector/sparse.h>

using namespace LELA;

template <class Ring>
void MyAlgorithm (const Ring &R)
{
	typename Vector<Ring>::Sparse v;

	// Index zero gets value -1
	v.push_back (typename Vector<Ring>::Sparse::value_type (0, Ring::Element ()));
	R.init (v.back ().second, -1);

	// Index 15 gets value 3
	v.push_back (typename Vector<Ring>::Sparse::value_type (15, Ring::Element ()));
	R.init (v.back ().second, 3);

	// Index 23 gets value 10
	v.push_back (typename Vector<Ring>::Sparse::value_type (23, Ring::Element ()));
	R.init (v.back ().second, 10);
}

#include <lela/ring/gf2.h>

using namespace LELA;

void MyAlgorithm (const GF2 &F, unsigned int n)
{
	Vector<GF2>::Dense v (n);

	*(v.word_begin ()) = 0xf0f0f0f0f0f0f0f0ULL;     // Set all of the first 64 bits of the vector simultaneously
	*(v.word_begin () + 1) = 0xaaaaaaaaaaaaaaaaULL; // Set all of the second 64 bits of the vector simultaneously
	v.back_word () = 0xeeee5555eeee5555ULL;         // Set the final 64 bits of the vector simultaneously
}

#include <lela/ring/gf2.h>

using namespace LELA;

void MyAlgorithm (const GF2 &F)
{
	Vector<GF2>::Sparse v;

	v.push_back (1);   // Index 1 gets value one
	v.push_back (15);  // Index 15 gets value one
	v.push_back (23);  // Index 23 gets value one
}

#include <lela/ring/gf2.h>

using namespace LELA;

void MyAlgorithm (const GF2 &R)
{
	typename Vector<GF2>::Hybrid v;

	// Append entries to the vector
	v.push_back (Vector<GF2>::Hybrid::value_type (0, 0xffff0000ffff0000ULL));  // First 64 indices
	v.push_back (Vector<GF2>::Hybrid::value_type (1, 0xf0f0f0f0f0f0f0f0ULL));  // Second 64 indices
	v.push_back (Vector<GF2>::Hybrid::value_type (5, 0xaaaa5555aaaa5555ULL));  // Indices 5 * 64 - 6 * 64
}

Matrices in LELA are pure containers. They know nothing about the underlying ring or its arithmetic. LELA provides two basic matrix-types: dense and sparse. The interfaces for these matrices are similar.

To declare a dense matrix, use DenseMatrix, which is parametrised by the element-type.

#include <lela/matrix/dense.h>

using namespace LELA;

template <class Ring>
void MyFunction (const Ring &R)
{
	DenseMatrix<typename Ring::Element> M (5, 5); // Declare a 5x5 matrix

	typename Ring::Element a;
	R.init (a, 5);
	M.setEntry (1, 2, a);  // Set entry (1,2) to a
}

Note that row- and column-indices start from zero.

Caution

The state of the entries of a newly defined dense matrix is undefined. In particular, all entries should be initialised before the dense matrix is used in any calculation. This can be done: with setEntry as above, using the row- or column-iterators as described below, by copying from an existing matrix with the BLAS-interface (see the section called “Vector- and matrix-arithmetic”), by scaling the matrix to zero with the BLAS-interface (see the section called “Vector- and matrix-arithmetic”), by reading the matrix from a file, or by initialising the matrix with a VectorStream (see reference-documentation).

Dense matrices provide three types of iterators: row, column, and raw. The row- and column-iterators iterate over the row- and column-vectors of the matrix, providing them as dense vectors.

#include <lela/matrix/dense.h>

using namespace LELA;

template <class Ring>
void MyFunction (const Ring &R)
{
	DenseMatrix<typename Ring::Element> M (5, 5); // Declare a 5x5 matrix

	// Iterate through all the rows, setting the entry 2 on each one to 5 in the ring
	typename DenseMatrix<typename Ring::Element>::RowIterator i;

	for (i = M.rowBegin (); i != M.rowEnd (); ++i)
		R.init ((*i)[2], 5);

	// Iterate through all the columns, setting the entry 1 on each one to 10 in the ring
	typename DenseMatrix<typename Ring::Element>::ColIterator j;

	for (j = M.colBegin (); j != M.colEnd (); ++j)
		R.init ((*i)[1], 10);
}

The raw iterator iterates through all the entries of the matrix in some order. It is useful for applying an operation uniformly to all entries in the matrix.

#include <lela/matrix/dense.h>

using namespace LELA;

template <class Ring>
void MyFunction (const Ring &R, DenseMatrix<typename Ring::Element> &M)
{
	// Set all entries in the matrix to 5 in the ring
	DenseMatrix<typename Ring::Element>::RawIterator i;

	for (i = M.rawBegin (); i != M.rawEnd (); ++i)
		R.init (*i, 5);
}

Each of the iterators has a const equivalent, the name prefixed by Const.

Over GF(2), the declaration of a dense matrix is analogous: use DenseMatrix<bool>. There is however no support for column-iterators.

	Note
	If LELA is configured with support for M4RI, then DenseMatrix<bool> is a wrapper for the matrix-type in M4RI.

The sparse matrix stores the rows of the matrix as sparse vectors. Its type is SparseMatrix and it is parametrised by element-type.

#include <lela/matrix/sparse.h>

using namespace LELA;

template <class Ring>
void MyFunction (const Ring &R)
{
	SparseMatrix<typename Ring::Element> M (5, 5); // Declare a 5x5 matrix

	typename Ring::Element a;
	R.init (a, 5);
	M.setEntry (1, 2, a);  // Set entry (1,2) to a
}

Sparse matrices have no support for column-iterators. Row iterators iterate over row-vectors which are of a sparse type. Otherwise the interface is identical to that of DenseMatrix.

In contrast to dense matrices, a newly defined sparse matrix is by definition the zero-matrix.

Sparse matrices are also parameterised by the type of the row-vector. This can be used to define a sparse matrix over GF(2) whose rows are hybrid vectors.

#include <lela/ring/gf2.h>
#include <lela/matrix/sparse.h>

using namespace LELA;

void MyFunction (const GF2 &F)
{
	SparseMatrix<bool, Vector<GF2>::Hybrid> M (5, 5); // Declare a 5x5 matrix
	...
}

LELA provides support for declaring subvectors, submatrices, and transposed matrices which are just windows into their respective containers and therefore require no additional copying. There are, however, some caveats which we will discuss here.

There are two types of subvectors and respectively submatrices: arbitrary and aligned. The starting and possibly ending indices of an aligned subvector or submatrix may be limited to an integral multiple of some value. However, calculations with aligned subvectors and submatrices may be much faster than with arbitrary ones. In addition, the hybrid vector-representation does not permit arbitrary subvectors which can be modified. Aligned subvector and submatrices are useful in recursive algorithms which divide matrices and vectors up and where the decision of where to make the division is open to the implementor of the algorithm. Currently these types are only distinct with vectors over GF(2).

One may declare an arbitrary subvector of a vector of type Vector by defining an object of the type VectorTraits<Ring, Vector>::SubvectorType and passing to its constructor the parent-vector and the starting and ending indices of the subvector. Here we define a subvector w of a previously defined vector v to span the indices 10 through 100 (not including 100 itself).

Vector v;
...
typename VectorTraits<Ring, Vector>::SubvectorType w (v, 10, 100);

If the parent-vector is const, then it is required to use VectorTraits<Ring, Vector>::ConstSubvectorType instead.

const Vector v;
...
typename VectorTraits<Ring, Vector>::ConstSubvectorType w (v, 10, 100);

To declare an aligned subvector, use VectorTraits<Ring, Vector>::AlignedSubvectorType, respectively VectorTraits<Ring, Vector>::ConstAlignedSubvectorType. In this case the starting index must be an integral multiple of the static integer VectorTraits<Ring, Vector>::align. The ending index should either be an integral multiple of VectorTraits<Ring, Vector>::align or coincide with the end of the vector.

This basic mechanism works for all vector-types but there are a few limitations:

Caution

It is important to be extremely careful with mutable subvectors of sparse or hybrid vectors. If a subvector is changed, then other subvectors of the same vector will not be updated to reflect that change. For example, suppose you have two subvectors x_1 and x_2 of a vector x. If you append an entry to x_1 it will appear correctly in x, but x_2 will be left in an invalid state: Vector<Ring>::Sparse x; ... typename Vector<Ring>::Sparse::SubvectorType x_1 (x, 0, 10); typename Vector<Ring>::Sparse::SubvectorType x_2 (x, 10, 20); x_1.push_back (5, typename Ring::Element ()); // x_2 now invalid Failure to take this into account may lead to random memory-corruption. This can also be a problem when using the arithmetic-routines described below. It is advised to avoid using two subvectors of a given sparse vector as inputs to the same operation if one of them will be changed, for example: BLAS2::gemv (ctx, a, A, x_2, b, x_1); // DANGER: Using both x_1 and x_2!

Declaring a submatrix is similar. To create an arbitrary mutable submatrix of a given matrix of type Matrix, use Matrix::SubmatrixType and supply to the constructor the containing matrix, the row- and column-indices of the upper left entry (indexed from 0, as usual), and the row- and column-dimensions. The following creates a 5x5 submatrix starting at index (2,2) of a 10x10 matrix.

Matrix M (10, 10);  // M has size 10x10

typename Matrix::SubmatrixType S (M, 2, 2, 5, 5);

Note the difference with subvectors: for submatrices one supplies the dimensions of the submatrix, while for subvectors one supplies instead the ending index.

If the matrix is const, then use Matrix::ConstSubmatrixType.

const Matrix M (10, 10);  // M has size 10x10

typename Matrix::ConstSubmatrixType S (M, 2, 2, 5, 5);

One can create aligned submatrices with Matrix::AlignedSubmatrixType and Matrix::ConstAlignedSubmatrixType. In this case the starting row- and column-indices must be integer multiples of Matrix::rowAlign and respectively Matrix::colAlign. The dimensions should be such that either the ending row and column are aligned with an integer multiple of the above values or they coincide with the end of the matrix.

As with subvectors, it is not permitted to take an mutable arbitrary submatrix of a sparse matrix whose rows or columns are hybrid vectors. Unlike subvectors it is okay to construct a mutable submatrix of a mutable submatrix of a sparse matrix.

To construct a transposed matrix, use the class TransposeMatrix, defined in lela/matrix/transpose.h. This class is parametrised by matrix-type.

#include "lela/matrix/transpose.h"
...
Matrix M;

TransposeMatrix<Matrix> MT (M);

If the matrix is const, then TransposeMatrix should also be instantiated with the const variant.

const Matrix M;

TransposeMatrix<const Matrix> MT (M);

	Note
	In LELA there also exist classes Subvector, SparseSubvector, and Submatrix. These are basically for internal use. Do not use them unless you are creating your own vector- or matrix-type and wish to employ them to define the subvector or submatrix-types used above.

To compute the row echelon form of a matrix, use the solution EchelonForm. Its header-file is lela/solutions/echelon-form.h.

The use is quite straightforward:

#include <lela/solutions/echelon-form.h>

using namespace LELA;

template <class Ring, class Matrix>
void ComputeRowEchelonForm (const Ring &R, Matrix &A)
{
	Context<Ring> ctx (R);
	EchelonForm<Ring> EF (ctx);

	echelonize (A); // A now replaced by its row echelon form
}

There is also support for computing the reduced row echelon form of a matrix:

#include <lela/solutions/echelon-form.h>

using namespace LELA;

template <class Ring, class Matrix>
void ComputeRowEchelonForm (const Ring &R, Matrix &A)
{
	Context<Ring> ctx (R);
	EchelonForm<Ring> EF (ctx);

	echelonize (A, true); // A now replaced by its reduced row echelon form
}

You can also instruct EchelonForm to use a particular method when computing the row-echelon form. For example, if your matrix comes from the F4 algorithm, it is useful to use the method of Faugère-Lachartre:

#include <lela/solutions/echelon-form.h>

using namespace LELA;

template <class Ring, class Matrix>
void ComputeRowEchelonForm (const Ring &R, Matrix &A)
{
	Context<Ring> ctx (R);
	EchelonForm<Ring> EF (ctx);

	echelonize (A, true, METHOD_FAUGERE_LACHARTRE); // A now replaced by its reduced row echelon form
}

The supported methods are as follows:

Throughout LELA there are checks of the form

lela_check (condition);

These check that the input to algorithms meets the required preconditions. These checks are disabled by default, but can be enabled by defining the macro DEBUG when compiling both the library and sources.

LELA has a facility called the commentator for reporting information about computations it performs. This is particularly useful for obtaining benchmarks.

A program-run is divided conceptually into activities which form a tree. To mark the beginning and end of an activity, use the interfaces commentator.start and respectively commentator.stop.

#include <lela/util/commentator.h>

using namespace LELA;

void MyFunction ()
{
	commentator.start ("Description of activity", __FUNCTION__);
	...
	commentator.stop (MSG_DONE);
}

The second argument to commentator.start should be the name of the calling function. The argument to commentator.stop is a short string indicating the state at the end of the activity, such as whether a test passed. The macros MSG_DONE and MSG_STATUS provide standardised strings for this purpose.

One can use the commentator to report information about a computation with the interface commentator.report, which returns an std::ostream.

#include <lela/util/commentator.h>

using namespace LELA;

void MyFunction ()
{
	commentator.start ("Description of activity", __FUNCTION__);

	std::ostream &report = commentator.report (Commentator::LEVEL_IMPORTANT, INTERNAL_DESCRIPTION);

	report << "Some information about the computation." << std::endl;

	commentator.stop (MSG_DONE);
}

The first parameter to commentator.report indicates the importance of the information; it can be either Commentator::LEVEL_IMPORTANT, Commentator::LEVEL_NORMAL, or Commentator::LEVEL_UNIMPORTANT. The second parameter indicates the general category of the information. Possible values include INTERNAL_DESCRIPTION, INTERNAL_WARNING, INTERNAL_ERROR, TIMING_MEASURE, and PROGRESS_REPORT.

The commentator includes a sophisticated system to select which messages are printed and which are not. Messages to be printed can be selected according to the nesting-depth of the activity, the importance-level, and the classification. See the commentator's documentation for more information.

The commentator can produce both a brief output for the console and a more detailed output which can be directed to a file. See the commentator's documentation for more information.

To disable the commentator entirely, define the macro DISABLE_COMMENTATOR. Then all commentator-lines compile to nothing.