0

素数を法とする方程式で作業するために、フロートでガウス消去法を実行する一連の線形方程式を解くJavaでコードを書き直そうとしています。問題は、それが機能しないことであり、何が問題なのかわかりません。方程式の小さなセットでは機能するようですが、大きなセットでは機能しないため、デバッグが困難になります。

私のアルゴリズムは最初の行を取り、最初の要素の逆数を見つけることでこれを正規化し、行のすべての要素にこの逆数を掛けます。次に、この行を他の行から十分な回数減算して、最初の要素をゼロにします。次の反復では、次の行に移動し、行 i のピボット要素が列 i になるまで同じ手順を実行します。最後に、前の行からすべての行を減算して、すべての列 (最後のものを除く) のゼロ以外の要素を 1 つだけ作成します。(今のところ、私は double を使用していますが、これは必要ありませんが、これは問題にはなりません)。これが私のコードです：

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

これは小さな例 (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}};

どちらが正しいか (最初の変数 = 0、2 番目の変数 = 1.0、3 番目の変数 = 0)、WolframAlpha で 0*k^0 + 1*k^1 + 0*k^2 for k = 1..3 を確認できます。

この例では、10 個の変数と方程式 a*k^0 + b*k^1 + c*k^2... (mod 29) for 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で確認できます）。正解は (abc ...) = (8 13 9 13 4 27 18 10 12 24 15) です。

誰かが私の間違いを見つけることができますか? または、Gauss mod p の実行方法を誤解していますか?

4

1 に答える 1

0

最初の段階で i 番目の位置にゼロ以外の行が見つからないようです。inverseそれと、あなたの機能がどのように機能するのか私にはわからないという事実を考えると、そこで問題が発生している場合とそうでない場合があります。

第 2 段階で「ピボット」を探している理由がわかりません。最初の段階から彼らがどこにいるか知っています。実際、なぜ第 2 段階があるのか​​まったくわかりません。代わりに、最初の段階ですべての除去を行います。これにより、コードが大幅に明確になります。

doubleマトリックスで s を使用している理由がわかりません。また、なぜあちこちで使用しているのかわかりませんMath.abs。ここでは等値比較が適切です。

最後に、Vandermonde システムを解いています。それが単なるテスト ケースではなく、アプリケーションである場合は、おそらく代わりにラグランジュ補間を使用する必要があります。

于 2013-11-17T16:09:56.743 に答える