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?
