
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?