# HG changeset patch
# User Chris Hall
# 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/sage/matrix/matrix2.pyx Sun May 30 02:49:06 2010 -0700
+++ b/sage/matrix/matrix2.pyx Thu May 27 18:39:27 2010 -0600
@@ -4529,6 +4529,140 @@
self.cache('echelon_form', self)
verbose('done with gauss echelon form', tm)
+ def weak_popov_form(self, ascend=True):
+ """
+ This function computes a weak Popov form of a matrix over a rational
+ function field `k(x)`, for `k` a field.
+
+ INPUT:
+
+ - `ascend` - if True, rows of output matrix `W` are sorted so
+ degree (= the maximum of the degrees of the elements in
+ the row) increases monotonically, and otherwise degrees decrease.
+
+ OUTPUT:
+
+ A 3-tuple `(W,N,d)` consisting of two matrices over `k(x)` and a list
+ of integers:
+
+ 1. `W` - matrix giving a weak the Popov form of self
+ 2. `N` - matrix representing row operations used to transform
+ `self` to `W`
+ 3. `d` - degree of respective columns of W; the degree of a column is
+ the maximum of the degree of its elements
+
+ `N` is invertible over `k(x)`. These matrices satisfy the relation
+ `N*self = W`.
+
+ EXAMPLES:
+
+ The routine expects matrices over the rational function field, but
+ other examples below show how one can provide matrices over the ring
+ of polynomials (whose quotient field is the rational function field).
+
+ ::
+
+ sage: R. = GF(3)['t']
+ sage: K = FractionField(R)
+ sage: M = matrix([[(t-1)^2/t],[(t-1)]])
+ sage: M.weak_popov_form()
+ (
+ [ 0] [ t 2*t + 1]
+ [(2*t + 1)/t], [ 1 2], [-Infinity, 0]
+ )
+
+ If `self` is an `n x 1` matrix with at least one non-zero entry, `W` has
+ a single non-zero entry and that entry is a scalar multiple of the
+ greatest-common-divisor of the entries of `self`.
+
+ ::
+
+ sage: M = matrix([[t*(t-1)*(t+1)],[t*(t-2)*(t+2)],[t]])
+ sage: M.weak_popov_form()
+ (
+ [0] [ 1 0 2*t^2 + 1]
+ [0] [ 0 1 2*t^2 + 1]
+ [t], [ 0 0 1], [-Infinity, -Infinity, 1]
+ )
+
+ The following is the first half of example 5 in [H] *except* that we
+ have transposed `self`; [H] uses column operations and we use row.
+
+ ::
+
+ sage: R. = QQ['t']
+ sage: M = matrix([[t^3 - t,t^2 - 2],[0,t]]).transpose()
+ sage: M.weak_popov_form()
+ (
+ [ t -t^2] [ 1 -t]
+ [t^2 - 2 t], [ 0 1], [2, 2]
+ )
+
+ The next example demonstrates what happens when `self` is a zero matrix.
+
+ ::
+
+ sage: R. = GF(5)['t']
+ sage: K = FractionField(R)
+ sage: M = matrix([[K(0),K(0)],[K(0),K(0)]])
+ sage: M.weak_popov_form()
+ (
+ [0 0] [1 0]
+ [0 0], [0 1], [-Infinity, -Infinity]
+ )
+
+ In the following example, `self` has more rows than columns.
+
+ ::
+
+ sage: R. = QQ['t']
+ sage: M = matrix([[t,t,t],[0,0,t]], ascend=False)
+ sage: M.weak_popov_form()
+ (
+ [t t t] [1 0]
+ [0 0 t], [0 1], [1, 1]
+ )
+
+ The next example shows that M must be a matrix with
+ coefficients in a rational function field `k(t)`.
+
+ ::
+
+ sage: M = matrix([[1,0],[1,1]])
+ sage: M.weak_popov_form()
+ Traceback (most recent call last):
+ ...
+ TypeError: the coefficients of M must lie in a univariate
+ polynomial ring
+
+
+ NOTES:
+
+ - For consistency with LLL and other algorithms in sage, we have opted
+ for row operations; however, references such as [H] transpose and use
+ column operations.
+
+ - There are multiple weak Popov forms of a matrix, so one may want to
+ extend this code to produce a Popov form (see section 1.2 of [V]). The
+ latter is canonical, but more work to produce.
+
+ REFERENCES:
+
+ .. [H] F. Hess, "Computing Riemann-Roch spaces in algebraic function
+ fields and related topics," J. Symbolic Comput. 33 (2002), no. 4,
+ 425--445.
+
+ .. [MS] T. Mulders, A. Storjohann, "On lattice reduction for polynomial
+ matrices," J. Symbolic Comput. 35 (2003), no. 4, 377--401
+
+ .. [V] G. Villard, "Computing Popov and Hermite forms of polynomial
+ matrices", ISSAC '96: Proceedings of the 1996 international symposium
+ on Symbolic and algebraic computation, 250--258.
+
+ """
+ import sage.matrix.matrix_misc
+ return sage.matrix.matrix_misc.weak_popov_form(self)
+
#####################################################################################
# Windowed Strassen Matrix Multiplication and Echelon
# Precise algorithms invented and implemented by David Harvey and Robert Bradshaw
diff -r 5e31911eb5c3 -r 42c8c73476a6 sage/matrix/matrix_misc.py
--- a/sage/matrix/matrix_misc.py Sun May 30 02:49:06 2010 -0700
+++ b/sage/matrix/matrix_misc.py Thu May 27 18:39:27 2010 -0600
@@ -21,3 +21,173 @@
def row_iterator(A):
for i in xrange(A.nrows()):
yield A.row(i)
+
+def weak_popov_form(M,ascend=True):
+ """
+ This function computes a weak Popov form of a matrix over a rational
+ function field `k(x)`, for `k` a field.
+
+ INPUT:
+
+ - `M` - matrix
+
+ - `ascend` - if True, rows of output matrix `W` are sorted so
+ degree (= the maximum of the degrees of the elements in
+ the row) increases monotonically, and otherwise degrees decrease.
+
+ OUTPUT:
+
+ A 3-tuple `(W,N,d)` consisting of two matrices over `k(x)` and a list
+ of integers:
+
+ 1. `W` - matrix giving a weak the Popov form of M
+ 2. `N` - matrix representing row operations used to transform
+ `M` to `W`
+ 3. `d` - degree of respective columns of W; the degree of a column is
+ the maximum of the degree of its elements
+
+ `N` is invertible over `k(x)`. These matrices satisfy the relation
+ `N*M = W`.
+
+ EXAMPLES:
+
+ The routine expects matrices over the rational function field, but
+ other examples below show how one can provide matrices over the ring
+ of polynomials (whose quotient field is the rational function field).
+
+ ::
+
+ sage: R. = GF(3)['t']
+ sage: K = FractionField(R)
+ sage: import sage.matrix.matrix_misc
+ sage: sage.matrix.matrix_misc.weak_popov_form(matrix([[(t-1)^2/t],[(t-1)]]))
+ (
+ [ 0] [ t 2*t + 1]
+ [(2*t + 1)/t], [ 1 2], [-Infinity, 0]
+ )
+
+ NOTES:
+
+ See docstring for weak_popov_form method of matrices for
+ more information.
+ """
+
+ # determine whether M has polynomial or rational function coefficients
+ R = M.base_ring()
+ from sage.rings.ring import is_Field
+ if is_Field(R):
+ R = R.base()
+ from sage.rings.polynomial.polynomial_ring import is_PolynomialRing
+ if not is_PolynomialRing(R):
+ raise TypeError("the coefficients of M must lie in a univariate polynomial ring")
+
+ t = R.gen()
+
+ # calculate least-common denominator of matrix entries. we coerce it
+ # to lie in R in the event that the entries of M lie in R (hence
+ # have no denominator)
+ from sage.rings.arith import lcm
+ den = R(lcm([a.denominator() for a in M.list()]))
+
+ # clear denominators
+ from sage.matrix.constructor import matrix
+ from sage.misc.functional import numerator
+ num = matrix([(lambda x : map(numerator, x))(v) for v in map(list,(M*den).rows())])
+
+ r = [list(v) for v in num.rows()]
+
+ N = matrix(num.nrows(), num.nrows(), R(1)).rows()
+
+ from sage.rings.infinity import Infinity
+ if M.is_zero():
+ return (M, matrix(N), [-Infinity for i in range(num.nrows())])
+
+ rank = 0
+ num_zero = 0
+ while rank != len(r) - num_zero:
+ # construct matrix of leading coefficients
+ v = []
+ for w in map(list, r):
+ # calculate degree of row (= max of degree of entries)
+ d = max([e.numerator().degree() for e in w])
+
+ # extract leading coefficients from current row
+ x = []
+ for y in w:
+ if y.degree() >= d and d >= 0: x.append(y.coeffs()[d])
+ else: x.append(0)
+ v.append(x)
+ l = matrix(v)
+
+ # count number of zero rows in leading coefficient matrix
+ # because they do *not* contribute interesting relations
+ num_zero = 0
+ for v in l.rows():
+ is_zero = 1
+ for w in v:
+ if w != 0:
+ is_zero = 0
+ if is_zero == 1:
+ num_zero += 1
+
+ # find non-trivial relations among the columns of the
+ # leading coefficient matrix
+ kern = l.kernel().basis()
+ rank = num.nrows() - len(kern)
+
+ # do a row operation if there's a non-trivial relation
+ if not rank == len(r) - num_zero:
+ for rel in kern:
+ # find the row of num involved in the relation and of
+ # maximal degree
+ indices = []
+ degrees = []
+ for i in range(len(rel)):
+ if rel[i] != 0:
+ indices.append(i)
+ degrees.append(max([e.degree() for e in r[i]]))
+
+ # find maximum degree among rows involved in relation
+ max_deg = max(degrees)
+
+ # check if relation involves non-zero rows
+ if max_deg != -1:
+ i = degrees.index(max_deg)
+ rel /= rel[indices[i]]
+
+ for j in range(len(indices)):
+ if j != i:
+ # do row operation and record it
+ v = []
+ for k in range(len(r[indices[i]])):
+ v.append(r[indices[i]][k] + rel[indices[j]] * t**(max_deg-degrees[j]) * r[indices[j]][k])
+ r[indices[i]] = v
+
+ v = []
+ for k in range(len(N[indices[i]])):
+ v.append(N[indices[i]][k] + rel[indices[j]] * t**(max_deg-degrees[j]) * N[indices[j]][k])
+ N[indices[i]] = v
+
+ # remaining relations (if any) are no longer valid,
+ # so continue onto next step of algorithm
+ break
+
+ # sort the rows in order of degree
+ d = []
+ from sage.rings.all import infinity
+ for i in range(len(r)):
+ d.append(max([e.degree() for e in r[i]]))
+ if d[i] < 0:
+ d[i] = -infinity
+ else:
+ d[i] -= den.degree()
+
+ for i in range(len(r)):
+ for j in range(i+1,len(r)):
+ if (ascend and d[i] > d[j]) or (not ascend and d[i] < d[j]):
+ (r[i], r[j]) = (r[j], r[i])
+ (d[i], d[j]) = (d[j], d[i])
+ (N[i], N[j]) = (N[j], N[i])
+
+ # return reduced matrix and operations matrix
+ return (matrix(r)/den, matrix(N), d)