Finite Field C++ Implementation

I am implementing finite field arithmetic for some research purposes in C++. The field of order v, when a prime (and not a prime power), is just modular arithmetic modulo v. Otherwise, v could be a prime power, where the arithmetic is not straightforward. For simplicity, assume that files that contain the multiplication and addition tables for all necessary prime powers are stored in the directory of this program. (I also need to compute the subtraction, division, additive inverse, and multiplicative inverse tables).

My strategy is to store the multiplication/addition tables as a 2D std::vector and when I want to compute i*j in the field, I just return the corresponding table lookup. I figured that since the tables will never change I could declare them const, but I need to guarantee an order on when the tables are computed initially, since computing some tables require others to be finished before they are computed.

The methods isFullRank, canExchange, eliminate, and eliminateRow are all performing a standard algorithm for determining if a given square matrix is singular (i.e., full rank, or determinant equals 0) or not over the field. I was not able to find any library that implements this over the field. Although the standard technique of determinant calculation and modulo the field order at the end works for prime v, it does not necessarily work for prime power v.

This is FiniteField.h:

#include <fstream>
#include <sstream>
#include <vector>

class FiniteField {
public:
FiniteField(int v) : v(v) {
this->mulTable = std::vector<std::vector<int>>(v, std::vector<int>(v, 0));
this->addTable = std::vector<std::vector<int>>(v, std::vector<int>(v, 0));
this->subTable = std::vector<std::vector<int>>(v, std::vector<int>(v, 0));
this->divTable = std::vector<std::vector<int>>(v, std::vector<int>(v, 0));
this->mulInverse = std::vector<int>(v);
this->computeMulTable();
this->computeMulInverseTable();
this->computeSubTable();
this->computeDivTable();
}

int sub(int, int);
int mul(int, int);
int div(int, int);

bool isFullRank(std::vector<std::vector<int>>&);
bool isPrime(int n);
private:
const int v;
std::vector<std::vector<int>> mulTable;
std::vector<std::vector<int>> subTable;
std::vector<std::vector<int>> divTable;
std::vector<int> mulInverse;

void computeMulTable();
void computeSubTable();
void computeDivTable();

void computeMulInverseTable();

int canExchange(std::vector<std::vector<int>>& a, const int i);
void eliminate(std::vector<std::vector<int>>& a, const int i);
void eliminateRow(std::vector<std::vector<int>>& a, const int i, const int j);
};

And here is FiniteField.cpp:

#include "FiniteField.h"

bool FiniteField::isPrime(int n) {
for (int i = 2; (i * i) <= n; i++) {
if (n % i == 0) {
return false;
}
}
return true;
}

void FiniteField::computeMulTable() {
if (isPrime(this->v)) {
for (int i = 0; i < this->v; i++) {
for (int j = 0; j < this->v; j++) {
this->mulTable[i][j] = (i * j) % this->v;
}
}
} else {
// load multiplication table from a file
std::string fileName = "mul" + std::to_string(this->v) + ".txt";
std::ifstream myFile(fileName);

if (myFile.good()) {
for (int i = 0; i < this->v; i++) {
std::vector<int> row;
std::string line;
std::getline(myFile, line);
auto ss = std::istringstream(line);
for (int j = 0; j < this->v; j++) {
int p;
ss >> p;
row.push_back(p);
}
this->mulTable[i] = row;
}
} else {
exit(1);
}
}
}

void FiniteField::computeMulInverseTable() {
for (int i = 0; i < this->v; i++) {
for (int j = 0; j < this->v; j++) {
if (this->mulTable[i][j] == 1) {
this->mulInverse[i] = j;
}
}
}
}

if (isPrime(this->v)) {
for (int i = 0; i < this->v; i++) {
for (int j = 0; j < this->v; j++) {
this->addTable[i][j] = (i + j) % this->v;
}
}
} else {
// load addition table from a file
std::string fileName = "add" + std::to_string(this->v) + ".txt";
std::ifstream myFile(fileName);

if (myFile.good()) {
for (int i = 0; i < this->v; i++) {
std::vector<int> row;
std::string line;
std::getline(myFile, line);
auto ss = std::istringstream(line);
for (int j = 0; j < this->v; j++) {
int p;
ss >> p;
row.push_back(p);
}
}
} else {
exit(1);
}
}
}

for (int i = 0; i < this->v; i++) {
for (int j = 0; j < this->v; j++) {
if (this->addTable[i][j] == 0) {
}
}
}
}

void FiniteField::computeSubTable() {
for (int i = 0; i < this->v; i++) {
for (int j = 0; j < this->v; j++) {
}
}
}

void FiniteField::computeDivTable() {
for (int i = 0; i < this->v; i++) {
for (int j = 0; j < this->v; j++) {
this->divTable[i][j] = this->mulTable[i][this->mulInverse[j]];
}
}
}

int FiniteField::add(int i, int j) {
}

int FiniteField::mul(int i, int j) {
return this->mulTable[i][j];
}

int FiniteField::sub(int i, int j) {
return this->subTable[i][j];
}

int FiniteField::div(int i, int j) {
return this->divTable[i][j];
}

bool FiniteField::isFullRank(std::vector<std::vector<int>>& a) {
int n = a.size();
for (int i = 0; i < n; i++) {
if (a[i][i] == 0) {
int j = this->canExchange(a, i);
if (j < i) {
return false;
} else {
std::swap(a[i], a[j]);
}
}
this->eliminate(a, i);
}
return true;
}

int FiniteField::canExchange(std::vector<std::vector<int>>&a, const int i) {
int n = a.size();
int fail = i - 1;
for (int j = i + 1; j < n; j++) {
if (a[j][i] != 0) {
return j;
}
}
return fail;
}

void FiniteField::eliminate(std::vector<std::vector<int>>&a, const int i) {
int n = a.size();
for (int j = i + 1; j < n; j++) {
if (a[j][i] != 0) {
this->eliminateRow(a, i, j);
}
}
}

void FiniteField::eliminateRow(std::vector<std::vector<int>>&a, const int i, const int j) {
int p = this->div(a[j][i], a[i][i]);
int n = a.size();
for (int k = i; k < n; k++) {
a[j][k] = this->sub(a[j][k], this->mul(p, a[i][k]));
}
}

I am mainly concerned about the performance of using this class in terms of the call to isFullRank (which indirectly relates to accessing the given addition and multiplication tables).

My compiler options are: g++ -std=c++11 -O3. I am able to use C++14 on this machine as well, so I'm open to suggestions on using some optimizations for that as well.

First, you should try to make your code const-correct. Any member function that doesn't modify *this should be qualified with const, and any member function that doesn't even access *this should be qualified with static:

int add(int, int) const;
int sub(int, int) const;
int mul(int, int) const;
int div(int, int) const;

bool isFullRank(std::vector<std::vector<int>>&) const;
static bool isPrime(int n);

Next, any pointer or reference whose pointee is never modified should be qualified with const:

bool isFullRank(const std::vector<std::vector<int>>&) const;

Obviously your isPrime isn't the most efficient, but it's not on the hot path, so whatever.

int FiniteField::canExchange(std::vector<std::vector<int>>&a, const int i)

First, avoid const-qualifying non-pointerish parameters if you can; all that does is prevent you from moving out of them in C++11. (Not that it matters for int; but it can matter for more complicated types, and meanwhile it's not buying you anything.) But do const-qualify the pointee of that first reference parameter.

int FiniteField::canExchange(const std::vector<std::vector<int>>& a, int i)

More importantly, what is this function doing? Its name suggests that it should be returning bool: either you canExchange or you can't, right? What's the point of returning an integer?

Reading the code of the function makes it clear that the function is definitely misnamed... but I can't figure out what its name is supposed to be.

Also, you return the magic value i-1 to signal failure, and then in the caller you test for j == i-1 with j < i. It would be much much simpler and easier to reason about if you returned 0 on failure and tested for j == 0. (It might also be just a hair faster.)

bool FiniteField::isFullRank(std::vector<std::vector<int>>& a)

I was going to say "const-qualify your reference parameter" here, too; but it looks like your function actually does modify a! This is a big no-no. You should make a copy of the matrix if you're going to modify it inside the function, so that the caller's original matrix doesn't suddenly contain different values just because he called isFullRank.

In C++11, the way you unconditionally make a copy of a parameter is to just drop the ampersand:

bool FiniteField::isFullRank(std::vector<std::vector<int>> a)

This way, if the caller actually doesn't care about the values in his original matrix anymore, he can move it into place for you:

if (isFullRank(std::move(myMatrix))) {  // myMatrix is dead from here on
...
}

If you're concerned about raw speed, consider whether you can nail down the exact value of v at compile time. Obviously the compiler can generate much better code for e.g. "multiplication by 17" than it can for "multiplication by this->v". If you can nail down v at compile time, then you can make the whole FiniteField class into a class template:

template<int V>
class FiniteField {
static int add(int, int);  // etc.
};

assert(FiniteField<17>::add(13, 14) == 10);

In C++14, you could even make a bunch of these functions constexpr, so that e.g.

static_assert(FiniteField<17>::add(13, 14) == 10);

I wouldn't recommend trying that with C++11's limited form of constexpr, though. And I wouldn't try it at all unless you were looking at the generated assembly code and seeing something that wasn't being optimized appropriately.

std::vector<std::vector<int>> is probably not the most efficient representation for an NxN matrix of integers. Did you consider creating your own matrix class template, e.g. my::matrix<int, 4, 4> fourbyfour; and then making these functions take my::matrix<int, N, N>?

Admittedly you're making good use of the vector-of-vectors representation by taking advantage of its cheap pointerwise row-swap operation. But depending on the size of these matrices, it might even be cheaper to swap the rows valuewise than to be following pointers on the heap all the time. (The only way to know is to code it both ways and measure.)

• Could std::array work as well for the "more efficient representation" you're talking about? – Ryan Dec 27 '16 at 2:30
• Yes, std::array<std::array<int, 4>, 4> would eliminate pointer-chasing, but its row-swap operation would be expensive. I'm not sure if there might be a middle-ground option. – Quuxplusone Dec 27 '16 at 2:35
• I was able to achieve a 20-30% performance improvement in the structure I'm studying that uses this class by using std::array among some of the other suggestions you gave. Thanks! – Ryan Dec 27 '16 at 4:18