## So, dear reader, I wrote a Gauss-Jordan solver

Been a while…

…Anyway: I am now designing (PhD related) a simulation of a NoC system and I need to have it run some code – so I thought I’d look at code to solve a system of linear equations and sought to build the same using Gauss-Jordan elimination.

It proved to be a lot more difficult than I had hoped: not that the maths was too complex, but that producing something I could later simplify from C++ into something that looked like synthetic assembly has and is proving very difficult indeed.

My system needs to execute a simple synthetic assembly because what I am modelling is the way the memory interactions work – so I have to get close to the metal and not hide behind a high-level language (no, not even C will do). Obviously – or at least it feels like that to me – I don’t want to have to mess about with floating point implementations, so I wanted integer arithmetic.

And that is where the first problem comes – to implement a Gauss-Jordan solver (or even to solve linear equations in the way you used to at school: the Gauss-Jordan algorithm is merely a scaling up of that), you have to divide numbers and sadly (for this purpose) not all numbers are multiples of one another – so the only way to preserve integer maths is to use rationals.

That, of course, is not such a problem. The C++ `pair` utility even makes it easy to model a rational out of two integers (writing a separate `rational` class would be overkill for something as small as this.)

But then, quickly, another problem becomes apparent. Even when using `long long` types to represent the numerator and denominator of our rational fraction repeated multiplications to enable fractional subtraction breaks the integer bounds after little more than a handful of row operations (and I want to model a large system – with literally thousands of equations in each system).

Even when adding a repeatedly applied `gcd` function to apply Euclid’s algorithm to the rational fractions, the numbers go out of range quickly.

So I have been forced to rely on an arbitrary precision library – GNU’s GMP. The code now works well (you can see at its Github repro and the solver code is below:

```#include <iostream>
#include <sstream>
#include <fstream>
#include <string>
#include <vector>
#include <cstdlib>
#include <utility>
#include <gmpxx.h>

//Gauss-Jordan elimination

using namespace std;

void gcd(pair<mpz_class, mpz_class>& input)
{
mpz_class first = abs(input.first);
mpz_class second = abs(input.second);
if (first > second) {
first = second;
second = abs(input.first);
}
while (second != 0) {
mpz_class temp = second;
second = first%second;
first = temp;
}
input.first /= first;
input.second /= first;
if (input.second < 0) {
input.first *= -1;
input.second *= -1;
}
}

int main()
{
string path("./variables.csv");
ifstream inputFile(path);

string number;
while(getline(stringy, number, ',')) {
}

vector<vector<pair<mpz_class, mpz_class> > > lines;
vector<pair<mpz_class, mpz_class> > innerLine;
while (getline(stringy, number, ',')) {
pair<mpz_class, mpz_class>
}
lines.push_back(innerLine);
}

for (int i = 0; i < lines.size(); i++) {
pair<mpz_class, mpz_class> pivot(lines[i][i].second,
lines[i][i].first);
if (lines[i][i].first != 0) {
lines[i][i].first = 1;
lines[i][i].second = 1;
} else {
continue;
}
for (int j = i + 1; j <= lines.size(); j++) {
lines[i][j].first *= pivot.first;
lines[i][j].second *= pivot.second;
gcd(lines[i][j]);
}
for (int j = i + 1; j < lines.size(); j++) {
pair<mpz_class, mpz_class> multiple(lines[j][i].first,
lines[j][i].second);
lines[j][i] = pair<mpz_class, mpz_class>(0, 1);
for (int k = i + 1; k <= lines.size(); k++) {
pair<mpz_class, mpz_class>
factor(
multiple.first * lines[i][k].first,
multiple.second * lines[i][k].second);
gcd(factor);
lines[j][k] = pair<mpz_class, mpz_class>
(lines[j][k].first * factor.second -
factor.first * lines[j][k].second,
lines[j][k].second * factor.second);
gcd(lines[j][k]);
}

}
}
// now eliminate upper half
for (int i = lines.size() - 1; i > 0; i--) {
if (lines[i][i].first == 0) {
continue;
}
for (int j = i - 1; j >= 0; j--) {
pair<mpz_class, mpz_class> multiple = lines[j][i];
if (multiple.first == 0) {
continue;
}
lines[j][i] = pair<mpz_class, mpz_class>(0, 1);
lines[j][lines.size()].first =
lines[j][lines.size()].first * multiple.second
lines[j][lines.size()].second;
lines[j][lines.size()].second =
lines[j][lines.size()].second *
gcd(lines[j][lines.size()]);
}
}
cout << "DIAGONAL FORM" << endl;
for (int i = 0; i < lines.size(); i++) {
for (int j = 0; j < lines.size(); j++) {
if (lines[i][j].first == 0) {
cout << "0 , ";
} else {
cout << lines[i][j].first << "/" << lines[i][j].second << " , ";
}
}
cout << " == " << lines[i][lines.size()].first << " / " << lines[i][lines.size()].second << endl;
}

}
```

(NB: I know that GMP supports rationals “out of the box” but I am trying to model something that I can translate to my synthetic assembly and so hiding the details of the implementation is not what I want to do.)

How on Earth, though, am I going to model arbitrary precision code in my “assembly” code?

### One response to “So, dear reader, I wrote a Gauss-Jordan solver”

1. prubin73

This suggestion might be orthogonal to reality, but … if you have linear systems with exclusively integer coefficients and integer variables (and some reason to believe an all-integer solution exists), perhaps you should skip the linear algebra approach and mimic a constraint solver? You might not be able to get the speed you need, but it would involve exclusively integer arithmetic.