Гауссово исключение по модулю p

Question

Я пытаюсь переписать код в Java, решая набор линейных уравнений, выполняющих гауссовское исключение на поплавках, для работы с уравнениями по модулю. Проблема в том, что она не работает, и я не могу понять, что случилось. Кажется, он работает на небольших наборах уравнений, но не на больших наборах, что затрудняет отладку.

Мой алгоритм принимает первую строку, нормализует это путем нахождения инверсного первого элемента и умножает каждый элемент в строке на этот обратный. Затем он вычитает эту строку из других строк достаточно времени, чтобы сделать их первый элемент равным нулю. На следующей итерации он переходит к следующей строке и выполняет ту же процедуру, пока элемент поворота строки i не окажется в столбце i. Наконец, он вычитает каждую строку из предыдущих строк, чтобы сделать только один ненулевой элемент каждого столбца (кроме последнего). (На данный момент я использую двойники, что необязательно, но это не должно быть проблемой). Вот мой код:

// 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; 
     } 
    } 

Это работает на небольших примерах (мод 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}}; 

Что является правильным (первая переменная = 0, второй = 1,0, третья = 0), так как можно проверить по WolframAlpha для 0 * k^0 + 1 * k^1 + 0 * k^2 при k = 1..3.

Для этого примера, имеющего 10 переменных, и уравнения a * k^0 + b * k^1 + c * k^2 ... (mod 29) при k = 1..11, у меня это матрица:

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 

Используя мой алгоритм я получаю ответ:

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 

Но это неправильно! (можно проверить с помощью WolframAlpha). Правильный ответ должен быть (a b c ...) = (8 13 9 13 4 27 18 10 12 24 15).

Может ли кто-нибудь определить мою ошибку? Или я неправильно понял, как делать Gauss mod p?

Answer 1

Кажется, что вы не находите строку с ненулевым значением в i-й позиции на первом этапе. Учитывая это, и тот факт, что я не знаю, как работает ваша функция inverse, вы можете столкнуться с проблемой или не столкнуться с ней.

Я не понимаю, почему вы ищете «опорные точки» на втором этапе. Вы знаете, где они находятся на первом этапе. На самом деле, я не понимаю, почему у вас второй этап. Выполните все исключения на первом этапе. Это значительно упростит ваш код.

Я не понимаю, почему вы используете double s в своей матрице. Я также не понимаю, почему вы используете Math.abs повсюду; сравнений сравнений здесь достаточно.

Наконец, вы решаете систему Вандермонда. Если это ваше приложение, а не только тестовый пример, вы, вероятно, должны использовать интерполяцию Лагранжа.