StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

POGaussian elimination modulo p
text
Body
copied!<p>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. </p> <p>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:</p> <pre><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; } } </code></pre> <p>This works on small examples (mod 29):</p> <pre><code>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}}; </code></pre> <p>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.</p> <p>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:</p> <pre><code>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 </code></pre> <p>Using my algorithm I get the answer: </p> <pre><code>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 </code></pre> <p>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). </p> <p>Can anybody spot my mistake? Or have I misunderstood how to do Gauss mod p?</p>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload