Ticket #9069: trac_9069.patch

File trac_9069.patch, 12.3 KB (added by cjh, 11 years ago)
  • sage/matrix/matrix2.pyx

    # HG changeset patch
    # User Chris Hall <chall14@uwyo.edu>
    # Date 1275007167 21600
    # Node ID 42c8c73476a6fad6b24a47e6eb91a74c1178d740
    # Parent  5e31911eb5c3df4c54ec446126a8b126bc42e5f3
    Trac 9069: implementation of weak Popov form for matrices over a rational function field k(t)
    
    diff -r 5e31911eb5c3 -r 42c8c73476a6 sage/matrix/matrix2.pyx
    a b  
    45294529        self.cache('echelon_form', self)
    45304530        verbose('done with gauss echelon form', tm)
    45314531
     4532    def weak_popov_form(self, ascend=True):
     4533        """
     4534        This function computes a weak Popov form of a matrix over a rational
     4535        function field `k(x)`, for `k` a field.
     4536       
     4537        INPUT:
     4538           
     4539         - `ascend` - if True, rows of output matrix `W` are sorted so
     4540           degree (= the maximum of the degrees of the elements in
     4541           the row) increases monotonically, and otherwise degrees decrease.
     4542         
     4543        OUTPUT:
     4544       
     4545        A 3-tuple `(W,N,d)` consisting of two matrices over `k(x)` and a list
     4546        of integers:
     4547           
     4548        1. `W` - matrix giving a weak the Popov form of self
     4549        2. `N` - matrix representing row operations used to transform
     4550           `self` to `W`
     4551        3. `d` - degree of respective columns of W; the degree of a column is
     4552           the maximum of the degree of its elements
     4553       
     4554        `N` is invertible over `k(x)`. These matrices satisfy the relation
     4555        `N*self = W`.
     4556       
     4557        EXAMPLES:
     4558           
     4559        The routine expects matrices over the rational function field, but
     4560        other examples below show how one can provide matrices over the ring
     4561        of polynomials (whose quotient field is the rational function field).
     4562
     4563        ::
     4564               
     4565            sage: R.<t> = GF(3)['t']
     4566            sage: K = FractionField(R)
     4567            sage: M = matrix([[(t-1)^2/t],[(t-1)]])
     4568            sage: M.weak_popov_form()
     4569            (
     4570            [          0]  [      t 2*t + 1]               
     4571            [(2*t + 1)/t], [      1       2], [-Infinity, 0]
     4572            )
     4573                       
     4574        If `self` is an `n x 1` matrix with at least one non-zero entry, `W` has
     4575        a single non-zero entry and that entry is a scalar multiple of the
     4576        greatest-common-divisor of the entries of `self`.
     4577       
     4578        ::
     4579
     4580            sage: M = matrix([[t*(t-1)*(t+1)],[t*(t-2)*(t+2)],[t]])
     4581            sage: M.weak_popov_form()
     4582            (
     4583            [0]  [        1         0 2*t^2 + 1]                           
     4584            [0]  [        0         1 2*t^2 + 1]                           
     4585            [t], [        0         0         1], [-Infinity, -Infinity, 1]
     4586            )
     4587           
     4588        The following is the first half of example 5 in [H] *except* that we
     4589        have transposed `self`; [H] uses column operations and we use row.
     4590       
     4591        ::
     4592           
     4593            sage: R.<t> = QQ['t']
     4594            sage: M = matrix([[t^3 - t,t^2 - 2],[0,t]]).transpose()
     4595            sage: M.weak_popov_form()
     4596            (
     4597            [      t    -t^2]  [ 1 -t]       
     4598            [t^2 - 2       t], [ 0  1], [2, 2]
     4599            )
     4600                                   
     4601        The next example demonstrates what happens when `self` is a zero matrix.
     4602       
     4603        ::
     4604
     4605            sage: R.<t> = GF(5)['t']
     4606            sage: K = FractionField(R)
     4607            sage: M = matrix([[K(0),K(0)],[K(0),K(0)]])
     4608            sage: M.weak_popov_form()
     4609            (
     4610            [0 0]  [1 0]                       
     4611            [0 0], [0 1], [-Infinity, -Infinity]
     4612            )
     4613
     4614        In the following example, `self` has more rows than columns.
     4615
     4616        ::
     4617
     4618            sage: R.<t> = QQ['t']
     4619            sage: M = matrix([[t,t,t],[0,0,t]], ascend=False)
     4620            sage: M.weak_popov_form()
     4621            (
     4622            [t t t]  [1 0]       
     4623            [0 0 t], [0 1], [1, 1]
     4624            )
     4625
     4626        The next example shows that M must be a matrix with
     4627        coefficients in a rational function field `k(t)`.
     4628
     4629        ::
     4630
     4631            sage: M = matrix([[1,0],[1,1]])
     4632            sage: M.weak_popov_form()
     4633            Traceback (most recent call last):
     4634            ...
     4635            TypeError: the coefficients of M must lie in a univariate
     4636            polynomial ring
     4637
     4638
     4639        NOTES:
     4640       
     4641         - For consistency with LLL and other algorithms in sage, we have opted
     4642           for row operations; however, references such as [H] transpose and use
     4643           column operations.
     4644           
     4645         - There are multiple weak Popov forms of a matrix, so one may want to
     4646           extend this code to produce a Popov form (see section 1.2 of [V]).  The
     4647           latter is canonical, but more work to produce.
     4648                                   
     4649        REFERENCES:
     4650           
     4651        .. [H] F. Hess, "Computing Riemann-Roch spaces in algebraic function
     4652          fields and related topics," J. Symbolic Comput. 33 (2002), no. 4,
     4653          425--445.
     4654         
     4655        .. [MS] T. Mulders, A. Storjohann, "On lattice reduction for polynomial
     4656          matrices," J. Symbolic Comput. 35 (2003), no. 4, 377--401
     4657
     4658        .. [V] G. Villard, "Computing Popov and Hermite forms of polynomial
     4659          matrices", ISSAC '96: Proceedings of the 1996 international symposium
     4660          on Symbolic and algebraic computation, 250--258.
     4661         
     4662        """
     4663        import sage.matrix.matrix_misc
     4664        return sage.matrix.matrix_misc.weak_popov_form(self)
     4665
    45324666    #####################################################################################
    45334667    # Windowed Strassen Matrix Multiplication and Echelon
    45344668    # Precise algorithms invented and implemented by David Harvey and Robert Bradshaw
  • sage/matrix/matrix_misc.py

    diff -r 5e31911eb5c3 -r 42c8c73476a6 sage/matrix/matrix_misc.py
    a b  
    2121def row_iterator(A):
    2222    for i in xrange(A.nrows()):
    2323        yield A.row(i)
     24
     25def weak_popov_form(M,ascend=True):
     26    """
     27    This function computes a weak Popov form of a matrix over a rational
     28    function field `k(x)`, for `k` a field.
     29   
     30    INPUT:
     31       
     32     - `M` - matrix
     33     
     34     - `ascend` - if True, rows of output matrix `W` are sorted so
     35       degree (= the maximum of the degrees of the elements in
     36       the row) increases monotonically, and otherwise degrees decrease.
     37     
     38    OUTPUT:
     39   
     40    A 3-tuple `(W,N,d)` consisting of two matrices over `k(x)` and a list
     41    of integers:
     42       
     43    1. `W` - matrix giving a weak the Popov form of M
     44    2. `N` - matrix representing row operations used to transform
     45       `M` to `W`
     46    3. `d` - degree of respective columns of W; the degree of a column is
     47       the maximum of the degree of its elements
     48   
     49    `N` is invertible over `k(x)`. These matrices satisfy the relation
     50    `N*M = W`.
     51   
     52    EXAMPLES:
     53       
     54    The routine expects matrices over the rational function field, but
     55    other examples below show how one can provide matrices over the ring
     56    of polynomials (whose quotient field is the rational function field).
     57
     58    ::
     59           
     60        sage: R.<t> = GF(3)['t']
     61        sage: K = FractionField(R)
     62        sage: import sage.matrix.matrix_misc
     63        sage: sage.matrix.matrix_misc.weak_popov_form(matrix([[(t-1)^2/t],[(t-1)]]))
     64        (
     65        [          0]  [      t 2*t + 1]               
     66        [(2*t + 1)/t], [      1       2], [-Infinity, 0]
     67        )
     68
     69    NOTES:
     70
     71    See docstring for weak_popov_form method of matrices for
     72    more information.
     73    """
     74   
     75    # determine whether M has polynomial or rational function coefficients
     76    R = M.base_ring()
     77    from sage.rings.ring import is_Field
     78    if is_Field(R):
     79        R = R.base()
     80    from sage.rings.polynomial.polynomial_ring import is_PolynomialRing
     81    if not is_PolynomialRing(R):
     82        raise TypeError("the coefficients of M must lie in a univariate polynomial ring")
     83
     84    t = R.gen()
     85
     86    # calculate least-common denominator of matrix entries.  we coerce it
     87    # to lie in R in the event that the entries of M lie in R (hence
     88    # have no denominator)
     89    from sage.rings.arith import lcm
     90    den = R(lcm([a.denominator() for a in M.list()]))
     91
     92    # clear denominators
     93    from sage.matrix.constructor import matrix
     94    from sage.misc.functional import numerator
     95    num = matrix([(lambda x : map(numerator,  x))(v) for v in map(list,(M*den).rows())])
     96
     97    r = [list(v) for v in num.rows()]
     98
     99    N = matrix(num.nrows(), num.nrows(), R(1)).rows()
     100
     101    from sage.rings.infinity import Infinity
     102    if M.is_zero():
     103        return (M, matrix(N), [-Infinity for i in range(num.nrows())])
     104
     105    rank = 0
     106    num_zero = 0           
     107    while rank != len(r) - num_zero:
     108        # construct matrix of leading coefficients
     109        v = []
     110        for w in map(list, r):
     111            # calculate degree of row (= max of degree of entries)
     112            d = max([e.numerator().degree() for e in w])
     113
     114            # extract leading coefficients from current row
     115            x = []
     116            for y in w:
     117                if y.degree() >= d and d >= 0:   x.append(y.coeffs()[d])
     118                else:                            x.append(0)
     119            v.append(x)
     120        l = matrix(v)
     121       
     122        # count number of zero rows in leading coefficient matrix
     123        # because they do *not* contribute interesting relations
     124        num_zero = 0
     125        for v in l.rows():
     126            is_zero = 1
     127            for w in v:
     128                if w != 0:
     129                    is_zero = 0
     130            if is_zero == 1:
     131                num_zero += 1
     132       
     133        # find non-trivial relations among the columns of the
     134        # leading coefficient matrix
     135        kern = l.kernel().basis()
     136        rank = num.nrows() - len(kern)
     137
     138        # do a row operation if there's a non-trivial relation       
     139        if not rank == len(r) - num_zero:
     140            for rel in kern:
     141                # find the row of num involved in the relation and of
     142                # maximal degree           
     143                indices = []
     144                degrees = []
     145                for i in range(len(rel)):
     146                    if rel[i] != 0:
     147                        indices.append(i)
     148                        degrees.append(max([e.degree() for e in r[i]]))
     149               
     150                # find maximum degree among rows involved in relation       
     151                max_deg = max(degrees)
     152               
     153                # check if relation involves non-zero rows
     154                if max_deg != -1:
     155                    i = degrees.index(max_deg)
     156                    rel /= rel[indices[i]]
     157           
     158                    for j in range(len(indices)):
     159                        if j != i:
     160                            # do row operation and record it
     161                            v = []
     162                            for k in range(len(r[indices[i]])):
     163                                v.append(r[indices[i]][k] + rel[indices[j]] * t**(max_deg-degrees[j]) * r[indices[j]][k])
     164                            r[indices[i]] = v
     165
     166                            v = []
     167                            for k in range(len(N[indices[i]])):
     168                                v.append(N[indices[i]][k] + rel[indices[j]] * t**(max_deg-degrees[j]) * N[indices[j]][k])
     169                            N[indices[i]] = v
     170                   
     171                    # remaining relations (if any) are no longer valid,
     172                    # so continue onto next step of algorithm
     173                    break
     174
     175    # sort the rows in order of degree
     176    d = []
     177    from sage.rings.all import infinity
     178    for i in range(len(r)):
     179        d.append(max([e.degree() for e in r[i]]))
     180        if d[i] < 0:
     181            d[i] = -infinity
     182        else:
     183            d[i] -= den.degree()
     184       
     185    for i in range(len(r)):
     186        for j in range(i+1,len(r)):
     187            if (ascend and d[i] > d[j]) or (not ascend and d[i] < d[j]):
     188                (r[i], r[j]) = (r[j], r[i])
     189                (d[i], d[j]) = (d[j], d[i])
     190                (N[i], N[j]) = (N[j], N[i])
     191                           
     192    # return reduced matrix and operations matrix
     193    return (matrix(r)/den, matrix(N), d)