Revision history [back]

As stated in the question, one can obtain a pretty good result using the echelon form:

def mykernel(m):
    r = m.base_ring()
    m = m.echelon_form('frac')
    res = []
    for i in range(m.nrows(), m.ncols()):
        v = m.column(i)
        vr = v.base_ring()
        v = v.list()
        v = v + [0]*(m.ncols()-m.nrows())
        v[i] = -1
        v = vector(vr, v)
        v = v*v.denominator()
        v = v.change_ring(r)
        v = v/gcd(flatten([p.coefficients() for p in v]))
        res.append(v)
    return res

The code above makes assumptions which will be valid in most but not in all cases. In particular, it assumes that the left square of the echelon form will be an identity matrix. This is not neccessarily true, so a proper implementation would have to check for the actual location of the ones. But the general concept remains valid: every column in the echelon form which is not a unit vector corresponds to an element of the kernel.

In the example case, this yields the following vector:

[(-y1^2 + 4*y1*y2 - 4*y2^2 - 2*y1*y3 + 4*y2*y3 - y3^2, -x1^2 + 4*x1*x2 -
4*x2^2 - 2*x1*x3 + 4*x2*x3 - x3^2, 2*x1*y1 - 4*x2*y1 + 2*x3*y1 - 4*x1*y2
+ 8*x2*y2 - 4*x3*y2 + 2*x1*y3 - 4*x2*y3 + 2*x3*y3, 2*x3*y1^2 -
4*x2*y1*y2 - 4*x3*y1*y2 + 4*x1*y2^2 + 4*x3*y2^2 - 2*x1*y1*y3 +
8*x2*y1*y3 - 2*x3*y1*y3 - 4*x1*y2*y3 - 4*x2*y2*y3 + 2*x1*y3^2, 4*x2^2*y1
- 2*x1*x3*y1 - 4*x2*x3*y1 + 2*x3^2*y1 - 4*x1*x2*y2 + 8*x1*x3*y2 -
4*x2*x3*y2 + 2*x1^2*y3 - 4*x1*x2*y3 + 4*x2^2*y3 - 2*x1*x3*y3, -x3^2*y1^2
+ 4*x2*x3*y1*y2 - 4*x1*x3*y2^2 - 4*x2^2*y1*y3 + 2*x1*x3*y1*y3 +
4*x1*x2*y2*y3 - x1^2*y3^2)]