Gaussian elimination modulo p

Question

I'm trying to rewrite a code in Java solving a set of linear equations doing Gaussian Elimination on floats, to work on equations modulo a prime. The problem is that it does not work, and I cannot figure out what's wrong. It seems to work on small sets of equations, but not on large sets, which makes debugging difficult.

My algorithm take the first row, normalize this by finding the inverse to first element, and multiplies every element in the row by this inverse. Then it subtracts this row from the other rows enough times to make their first element zero. On the next iteration it goes to the next row and do the same procedure, until the pivoting element of row i is in column i. Lastly it subtracts every row from the previous rows to make only one non-zero element of every column (except the last one). (As of now I use doubles, which is not necessary, but this should not be a problem). Here's my code:

// Transforms A so that the leftmost square matrix has at most one 1 per row,
    // and no other nonzero elements.
    // O(n^3)
 public static void gauss(int[][] A, int num_columns) {
        int n = A.length;
        int m = A[0].length;

        for (int i = 0; i < num_columns; i++) {
            // Finding row with nonzero element at column i, swap this to row i
            for(int k = i; k < num_columns; k++){
                if(A[k][i] != 0){
                    int t[] = A[i];
                    A[i] = A[k];
                    A[k] = t;
                }
            }
            // Normalize the i-th row.
            int inverse = (int)inverse((long)A[i][i], prime);
            for (int k = i ; k < m; k++) A[i][k] = (A[i][k]*inverse) % prime;

            // Combine the i-th row with the following rows.
            for (int j = 0; j < n; j++) {
                if(j == i) continue;
                int c = A[j][i];
                A[j][i] = 0;
                for (int k = i + 1; k < m; k++){
                    A[j][k] = (A[j][k] - c * A[i][k] + c * prime) % prime;
                }
            }
        }
    }

    public static void gauss(int[][] A) {
        gauss(A, Math.min(A.length, A[0].length));
    }
    public static long gcd(long a, long b){
        if(a < b){
            long temp = a;
            a = b;
            b = temp;
        }
        if(b == 0) return a;
        return gcd(b, a % b);
     }
    public static Pair ext_euclid(long a, long b){
        if(a < b){
            Pair r = ext_euclid(b,a);
            return new Pair(r.second, r.first);
        }
        if(b == 0) return new Pair(1, 0);
        long q = a / b;
        long rem = a - b * q;
        Pair r = ext_euclid(b, rem);
        Pair ret = new Pair(r.second, r.first - q * r.second);
        return ret;
    }

    public static  long inverse(long num, long modulo){
        num = num%modulo;
        Pair p = ext_euclid(num, modulo);
        long ret = p.first;
        if(ret < 0) return (modulo + ret) % modulo;
        return ret % modulo;
    }

    static class Pair{
        public long first;
        public long second;
        public Pair(long frst, long scnd){
            first = frst;
            second = scnd;
        }
    }

This works on small examples (mod 29):

matrix = {{1.0, 1.0, 1.0, 1.0}, {1.0, 2.0, 1.0, 2.0},{1.0, 0.0, 0.0‚ 3.0}};
answer= {{1.0, 0.0, 0.0, 0.0},{0.0, 1.0, 0.0, 1.0}, {0.0, 0.0, 1.0, 0.0}};

Which is correct (first variable = 0, second = 1.0, third = 0), as can be checked by WolframAlpha for 0*k^0 + 1*k^1 + 0*k^2 for k = 1..3.

For this example, having 10 variables, and the equations a*k^0 + b*k^1 + c*k^2... (mod 29) for k = 1..11, I have this matrix:

1 1 1 1 1 1 1 1 1 1 1 8 
1 2 4 8 16 3 6 12 24 19 9 5 
1 3 9 27 23 11 4 12 7 21 5 12 
1 4 16 6 24 9 7 28 25 13 23 12 
1 5 25 9 16 22 23 28 24 4 20 15 
1 6 7 13 20 4 24 28 23 22 16 0 
1 7 20 24 23 16 25 1 7 20 24 5 
1 8 6 19 7 27 13 17 20 15 4 1 
1 9 23 4 7 5 16 28 20 6 7 18 
1 10 13 14 24 8 22 17 25 18 7 20 
1 11 5 26 25 14 9 12 16 7 7 8  

Using my algorithm I get the answer:

1 0 0 0 0 0 0 0 0 0 0 18 
0 1 0 0 0 0 0 0 0 0 0 8 
0 0 1 0 0 0 0 0 0 0 0 15 
0 0 0 1 0 0 0 0 0 0 0 11 
0 0 0 0 1 0 0 0 0 0 0 28 
0 0 0 0 0 1 0 0 0 0 0 27 
0 0 0 0 0 0 1 0 0 0 0 7 
0 0 0 0 0 0 0 1 0 0 0 21 
0 0 0 0 0 0 0 0 1 0 0 9 
0 0 0 0 0 0 0 0 0 1 0 24 
0 0 0 0 0 0 0 0 0 0 1 14  

But this is wrong! (can be checked with WolframAlpha). The right answer should be (a b c ...) = (8 13 9 13 4 27 18 10 12 24 15).

Can anybody spot my mistake? Or have I misunderstood how to do Gauss mod p?

Answer 1

You don't seem to find a row with a nonzero in the i'th position during your first stage. Given that, and the fact that I don't know how your inverse function works, you may or may not be running into a problem there.

I don't understand why you're looking for "pivots" in the second stage. You know where they are from the first stage. Actually, I don't see why you have a second stage at all. Do all of the elimination in the first stage instead. This will clarify your code greatly.

I don't understand why you're using doubles in your matrix. I also don't understand why you're using Math.abs all over the place; equality comparisons are plenty appropriate here.

Finally, you're solving a Vandermonde system. If that's your application, not just a test case, you should probably use Lagrange interpolation instead.