// compile and run with:
// c++ bareiss.cpp -lginac -lcln -ldl && ./a.out

#include <ginac/ginac.h>
#include <assert.h>

using namespace GiNaC;
using namespace std;

/* Since the internals of matrix class are protected, we need
 * to employ this hack to access them.
 */
class matrix_hack : public matrix {
    public:
    exvector &get_m() { return m; };
};

void
echelon_form_bareiss_original(matrix &m)
{
    unsigned nr = m.rows();
    unsigned nc = m.cols();
    exvector &mv = ((matrix_hack*)&m)->get_m();
    ex z = 1;
    unsigned r0 = 0;
    for (unsigned c0 = 0; (c0 < nc) && (r0 < nr - 1); c0++) {
        unsigned pivot = r0;
        // No normalization before is_zero() here, because
        // we maintain the matrix normalized throughout the
        // algorithm.
        while ((pivot < nr) && mv[pivot*nc + c0].is_zero()) pivot++;
        if (pivot == nr) {
            // The whole column below r0:c0 is zero, let's skip
            // it.
            continue;
        }
        if (pivot > r0) {
            // Found a non-zero row somewhere below r0; let's
            // swap it in.
            for (unsigned c = c0; c < nc; c++) {
                mv[pivot*nc + c].swap(mv[r0*nc + c]);
            }
        }
        ex a = mv[r0*nc + c0];
        for (unsigned r = r0 + 1; r < nr; r++) {
            ex b = mv[r*nc + c0];
            for (unsigned c = c0 + 1; c < nc; c++) {
                ex v;
                bool divok = divide(a*mv[r*nc + c] - b*mv[r0*nc + c], z, v);
                assert(divok);
                // Since v is a polynomial, the call to normal()
                // here does not involve GCD, but rather just
                // expands sums while preserving common factors,
                // which is nomally preferable to full expand().
                m[r*nc + c] = normal(v);
            }
            mv[r*nc + c0] = 0;
        }
        z = a;
        r0++;
    }
    // Zero out the remaining rows (just in case).
    for (unsigned r = r0 + 1; r < nr; r++) {
        for (unsigned c = 0; c < nc; c++) {
            mv[r*nc + c] = 0;
        }
    }
}

void
echelon_form_bareiss_skip_zrows(matrix &m)
{
    unsigned nr = m.rows();
    unsigned nc = m.cols();
    exvector &mv = ((matrix_hack*)&m)->get_m();
    ex z = 1;
    unsigned r0 = 0;
    for (unsigned c0 = 0; (c0 < nc) && (r0 < nr - 1); c0++) {
        unsigned pivot = r0;
        // No normalization before is_zero() here, because
        // we maintain the matrix normalized throughout the
        // algorithm.
        while ((pivot < nr) && mv[pivot*nc + c0].is_zero()) pivot++;
        if (pivot == nr) {
            // The whole column below r0:c0 is zero, let's skip
            // it.
            continue;
        }
        if (pivot > r0) {
            // Found a non-zero row somewhere below r0; let's
            // swap it in.
            for (unsigned c = c0; c < nc; c++) {
                mv[pivot*nc + c].swap(mv[r0*nc + c]);
            }
        }
        bool zerocol = true;
        for (unsigned r = r0 + 1; r < nr; r++) {
            ex b = mv[r*nc + c0];
            zerocol = zerocol && b.is_zero();
        }
        if (zerocol) {
            z = 1;
            r0++;
            continue;
        }
        ex a = mv[r0*nc + c0];
        for (unsigned r = r0 + 1; r < nr; r++) {
            ex b = mv[r*nc + c0];
            for (unsigned c = c0 + 1; c < nc; c++) {
                ex v;
                bool divok = divide(a*mv[r*nc + c] - b*mv[r0*nc + c], z, v);
                assert(divok);
                // Since v is a polynomial, the call to normal()
                // here does not involve GCD, but rather just
                // expands sums while preserving common factors,
                // which is nomally preferable to full expand().
                m[r*nc + c] = normal(v);
            }
            mv[r*nc + c0] = 0;
        }
        z = a;
        r0++;
    }
    // Zero out the remaining rows (just in case).
    for (unsigned r = r0 + 1; r < nr; r++) {
        for (unsigned c = 0; c < nc; c++) {
            mv[r*nc + c] = 0;
        }
    }
}

void
echelon_form_bareiss_skip_statfactors(matrix &m)
{
    unsigned nr = m.rows();
    unsigned nc = m.cols();
    exvector &mv = ((matrix_hack*)&m)->get_m();
    exvector extraf(nr);
    for (unsigned r = 0; r < nr; r++) { extraf[r] = 1; }
    ex z = 1;
    unsigned r0 = 0;
    for (unsigned c0 = 0; (c0 < nc) && (r0 < nr - 1); c0++) {
        unsigned pivot = r0;
        // No normalization before is_zero() here, because
        // we maintain the matrix normalized throughout the
        // algorithm.
        while ((pivot < nr) && mv[pivot*nc + c0].is_zero()) { pivot++; }
        if (pivot == nr) {
            // The whole column below r0:c0 is zero, we can skip
            // it.
            continue;
        }
        if (pivot > r0) {
            // Found a non-zero row somewhere below r0; let's
            // swap it in.
            for (unsigned c = c0; c < nc; c++) {
                mv[pivot*nc + c].swap(mv[r0*nc + c]);
            }
            extraf[pivot].swap(extraf[r0]);
        }
        ex a = mv[r0*nc + c0];
        for (unsigned r = r0 + 1; r < nr; r++) {
            ex b = mv[r*nc + c0];
            if (b.is_zero()) {
                // This branch is the same as the next one, but
                // is is simplified due to b being 0.
                ex f;
                assert(divide(a*extraf[r0]*extraf[r], z, f));
                extraf[r] = f;
            } else {
                ex g = gcd(a*extraf[r0], b*extraf[r]);
                ex f;
                gcd(extraf[r0]*g, z, &f);
                f = normal(f);
                ex d;
                assert(divide(z*f, extraf[r0]*extraf[r], d));
                for (unsigned c = c0 + 1; c < nc; c++) {
                    ex v;
                    assert(divide(a*mv[r*nc + c] - b*mv[r0*nc + c], d, v));
                    // Since v is a polynomial, the call to normal()
                    // here does not involve GCD, but rather just
                    // expands sums while preserving common factors,
                    // which is nomally preferable to full expand().
                    m[r*nc + c] = normal(v);
                }
                mv[r*nc + c0] = 0;
                extraf[r] = f;
            }
        }
        z = a*extraf[r0];
        r0++;
    }
    // Zero out the remaining rows (just in case).
    for (unsigned r = r0 + 1; r < nr; r++) {
        for (unsigned c = 0; c < nc; c++) {
            mv[r*nc + c] = 0;
        }
    }
    // At this point mv[r*nc + c]*extraf[r] is precisely what
    // plain Bareiss algorithm would have computed; we can now
    // discard this extra factor freely.
}

void
save_matrix(ostream &f, const matrix &m)
{
    f << "{";
    for (unsigned i = 0; i < m.rows(); i++) {
        if (i != 0) f << ",\n";
        f << "{";
        for (unsigned j = 0; j < m.cols(); j++) {
            if (j != 0) f << ", ";
            f << m(i, j);
        }
        f << "}";
    }
    f << "}";
}

int
main()
{
    symbol a11("a11"), a12("a12"), a13("a13"), a14("a14"), a15("a15");
    symbol a21("a21"), a22("a22"), a23("a23"), a24("a24"), a25("a25");
    symbol a31("a31"), a32("a32"), a33("a33"), a34("a34"), a35("a35");
    symbol a41("a41"), a42("a42"), a43("a43"), a44("a44"), a45("a45");
    symbol a51("a51"), a52("a52"), a53("a53"), a54("a54"), a55("a55");
    matrix m = {
        {a11, a12, a13, a14, a15},
        {a21, a22, a23, a24, a25},
        {  0, a32, a33, a34, a35},
        {  0,   0,   0, a44, a45},
        {  0,   0,   0,   0, a55}
    };
    cout << "original matrix:" << endl;
    save_matrix(cout, ex_to<matrix>(factor(m)));
    cout << endl << endl;

    matrix m1 = m;
    echelon_form_bareiss_original(m1);
    cout << "bareiss original:" << endl;
    save_matrix(cout, ex_to<matrix>(factor(m1)));
    cout << endl << endl;

    matrix m2 = m;
    echelon_form_bareiss_skip_zrows(m2);
    cout << "bareiss/skip zrows:" << endl;
    save_matrix(cout, ex_to<matrix>(factor(m2)));
    cout << endl << endl;

    matrix m3 = m;
    echelon_form_bareiss_skip_statfactors(m3);
    cout << "bareiss/skip statfactors:" << endl;
    save_matrix(cout, ex_to<matrix>(factor(m3)));
    cout << endl;
}