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); vector<mpz_class> answers; //answer line string rawAnswer; getline(inputFile, rawAnswer); istringstream stringy(rawAnswer); string number; while(getline(stringy, number, ',')) { answers.push_back(atol(number.c_str())); } //now read in the system vector<vector<pair<mpz_class, mpz_class> > > lines; while(getline(inputFile, rawAnswer)) { istringstream stringy(rawAnswer); vector<pair<mpz_class, mpz_class> > innerLine; while (getline(stringy, number, ',')) { pair<mpz_class, mpz_class> addPair(atol(number.c_str()), 1); innerLine.push_back(addPair); } 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; } pair<mpz_class, mpz_class> answer = lines[i][lines.size()]; 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 - answer.first * multiple.first * lines[j][lines.size()].second; lines[j][lines.size()].second = lines[j][lines.size()].second * answer.second * multiple.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”

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.