Ticket #174: trac-174.patch

File trac-174.patch, 6.2 KB (added by was, 14 years ago)
  • sage/matrix/matrix_integer_dense.pyx

    # HG changeset patch
    # User William Stein <wstein@gmail.com>
    # Date 1202450875 28800
    # Node ID 1662f89e140e7821b81bc5e26216bde00d1acb98
    # Parent  b28009b5acbd2e0524d6efc6d7020a3f42770acc
    trac #174 -- modular hermite normal form
    
    diff -r b28009b5acbd -r 1662f89e140e sage/matrix/matrix_integer_dense.pyx
    a b cdef class Matrix_integer_dense(matrix_d 
    24332433            return decomp_seq(D), decomp_seq(E)
    24342434        else:
    24352435            return decomp_seq([(W.intersection(V), t) for W, t in X])
     2436
     2437    def _add_row_and_maintain_echelon_form(self, row, pivots):
     2438        """
     2439        Assuming self is a full rank n x m matrix in reduced row
     2440        Echelon form over ZZ and row is a vector of degree m, this
     2441        function creates a new matrix that is the echelon form of self
     2442        with row appended to the bottom.
     2443
     2444        WARNING: It is assumed that self is in
     2445
     2446        INPUT:
     2447            row -- a vector of degree m over ZZ
     2448            pivots -- a list of integers that are the pivot columns of self.
     2449
     2450        OUTPUT:
     2451            matrix -- a matrix of in reduced row echelon form over ZZ
     2452            pivots -- list of integers
     2453
     2454        ALGORITHM: For each pivot column of self, we use the extended
     2455        Euclidean algorithm to clear the column.  The result is a new
     2456        matrix B whose row span is the same as self.stack(row), and
     2457        whose last row is 0 if and only if row is in the QQ-span of
     2458        the rows of self.  If row is not in the QQ-span of the rows of
     2459        self, then row is nonzero and suitable to be inserted into the
     2460        top n rows of A to form a new matrix that is in reduced row
     2461        echelon form.  We then clear that corresponding new pivot column.
     2462        """
     2463        cdef Py_ssize_t i, j, piv, n = self._nrows, m = self._ncols
     2464
     2465        import constructor
     2466
     2467        # 0. Base case
     2468        if self.nrows() == 0:
     2469            pos = row.nonzero_positions()
     2470            if len(pos) > 0:
     2471                pivots = [0]
     2472                i = pos[0]
     2473                if row[i] < 0:
     2474                    row *= -1
     2475            else:
     2476                pivots = []
     2477            return constructor.matrix([row]), pivots
     2478
    24362479       
     2480        # 1. Create a new matrix that has row as the last row.
     2481        row_mat = constructor.matrix(row)
     2482        A = self.stack(row_mat)
     2483       
     2484        # 2. Working from the left, clear each column to put
     2485        #    the resulting matrix back in echelon form.
     2486        for i, p in enumerate(pivots):
     2487            # p is the i-th pivot
     2488           
     2489            # (a). Take xgcd of pivot positions in last row and in ith
     2490            # row.
     2491
     2492            # TODO (optimize) -- change to use direct call to gmp and
     2493            # no bounds checking!
     2494            a = A[i,p]
     2495            b = A[n,p]
     2496            if b % a == 0:
     2497                # (b) Subtract a multiple of row i from row n.
     2498                c = b // a
     2499                if c:
     2500                    for j in range(m):
     2501                        A[n,j] -= c * A[i,j]
     2502            else:
     2503                # (b). More elaborate.
     2504                #  Replace the ith row by s*A[i] + t*A[n], which will
     2505                # have g in the i,p position, and replace the last row by
     2506                # (b//g)*A[i] - (a//g)*A[n], which will have 0 in the i,p
     2507                # position.
     2508                g, s, t = a.xgcd(b)
     2509
     2510                # TODO: change to use a single mpz_t* instead of an extracted row.
     2511                row_i = A.row(i)
     2512                row_n = A.row(n)
     2513
     2514                # TODO: change to use direct mpz_t access instead of []:
     2515                ag = a//g; bg = b//g
     2516
     2517                # OK -- now we have to make sure the top part of the matrix
     2518                # but with row i replaced by
     2519                #     r = s*row_i[j]  +  t*row_n[j]
     2520                # is put in rref.  We do this by recursively calling this
     2521                # function with the top part of A (all but last row) and the
     2522                # row r.
     2523
     2524                rows = A.rows()
     2525                new_top = s*row_i  +  t*row_n
     2526                new_bot = bg*row_i - ag*row_n
     2527               
     2528                del rows[-1]
     2529                del rows[i]
     2530                top_mat = constructor.matrix(ZZ, n-1, m, rows)
     2531                new_pivots = list(pivots)
     2532                del new_pivots[i]
     2533
     2534                top_mat, pivots = top_mat._add_row_and_maintain_echelon_form(new_top, new_pivots)
     2535                return top_mat._add_row_and_maintain_echelon_form(new_bot, pivots)
     2536
     2537        # 3. Insert last row in A sliding other rows down if it turns out
     2538        #     that the last row is nonzero.
     2539        # TODO: change to use some fast mpz_t access.
     2540        v = A.row(n)
     2541        new_pivots = list(pivots)
     2542        if v != 0:
     2543            R = A.rows()
     2544            # Determine where the last row should be inserted.
     2545            i = v.nonzero_positions()[0]
     2546            if i in pivots:
     2547                assert False, 'WARNING: bug in add_row -- i (=%s) should not be a pivot'%i
     2548            # If pivot entry is negative negate this row.
     2549            if v[i] < 0:
     2550                A.rescale_row(n, -1)
     2551                v = A.row(n)
     2552            new_pivots.append(i)
     2553            new_pivots.sort()
     2554            import bisect
     2555            j = bisect.bisect(pivots, i)
     2556            # The new row should go *before* row j, so it becomes row j
     2557            del R[-1]
     2558            R.insert(j, v)
     2559            A = A.parent()(R)
     2560
     2561        _clear_columns(A, new_pivots, A.nrows())
     2562
     2563        # end if
     2564        return A, new_pivots
     2565 
     2566def _clear_columns(A, pivots, n):
     2567    # Clear all columns
     2568    m = A.ncols()
     2569    for i, p in enumerate(pivots):
     2570        v = A.row(i)
     2571        b = v[p]
     2572        for k in range(n):
     2573            if k != i:
     2574                a = A[k,p]
     2575                if a:
     2576                    c = a//b
     2577                    # subtract off c*v from row k; resulting A[k,i] entry will be < b, hence in Echelon form.
     2578                    for l in range(m):
     2579                        A[k,l] -= c*v[l]
     2580                       
    24372581
    24382582###############################################################
    24392583