Lecture 42: Computing with Polynomials in Singular

Singular is a computer algebra system for polynomial computations, with special emphasis on commutative and noncommutative algebra, algebraic geometry, and singularity theory. Singular is integrated in Sage.

There are three ways to use Singular in Sage:

1. We can launch a Singular terminal window, or run a worksheet in Singular mode, so every command will be interpreted by Singular.

2. With singular() we have a Pythonic interface to Singular.

3. The input to singular.eval() is a string with a Singular command. The output is a string representation of a Singular object.

When we define a polynomial ring in Sage, then we are working with Singular objects.

Polynomials

We used Singular when we made polynomial rings.

R.<x,y> = PolynomialRing(QQ, order='lex')
R

In Sage we see Multivariate Polynomial Ring in x, y over Rational Field but a Singular user would get more technical information about the ring. We can get this as follows.

R._singular_()

and then we see as output

//   characteristic : 0
//   number of vars : 2
//        block   1 : ordering lp
//                  : names    x y
//        block   2 : ordering C

Polynomials in x and y automatically belong to this ring.

p = 4*y^2 - 4 + 4*x*y + x^2
q = 4*x^2 - 4*x*y + y^2 - 4
print(p, 'has type', type(p))
print(q, 'has type', type(q))

The type of both p and q is sage.rings.polynomial.multi_polynomial_libsingular.MPolynomial_libsingular. To see what we are solving for, we can plot the two curves defined by p and q.

cp = implicit_plot(p,(x,-2,2),(y,-2,2),color='red')
cq = implicit_plot(q,(x,-2,2),(y,-2,2),color='green')
(cp+cq).show(figsize=4)

The plot is shown in Fig. 97.

Fig. 97 Intersecting the green curve with the red curve.

We see that the curves break up in lines. Algebraically, we see the linear factors. Printing the outcome of factor(p) and factor(q) gives

(x + 2*y - 2) * (x + 2*y + 2)
(-2*x + y - 2) * (-2*x + y + 2)

Note that the objects we created actually belong to the Sage rings. We can declare Singular objects as follows.

sp = singular.new('4*y^2 - 4 + 4*x*y + x^2')
print(sp, 'has type', type(sp))

The type of sp is sage.interfaces.singular.SingularElement. Then on sp we can apply the sage() method

ssp = sp.sage()
print(ssp, 'has type', type(ssp))

The type of ssp is the same as p.

We can load a Singular library to solve polynomial systems.

singular.lib('solve.lib')
command = 'solve([%s, %s])' % (p, q)
print(command)
sols = singular.eval(command)
sols

From the output we see that we obtained four solutions. The four solutions are returned in a string, as the string representation of a list in Singular. We see that the solve in Singular solves over the complex numbers. What is returned is a string, as is always the case with the output of singular.eval().

Reducing Polynomials with a Groebner basis

The polynomials in the ring R define an ideal.

J = Ideal([p, q])
J

The ideal J consists of all polynomial combinations we can make with its generators p and q. And then we can compute a Groebner basis.

G = singular.groebner(J)
G

Because the term order was lexicographic, the result is triangular. We can solve x in terms of y.

singular.reduce(x,G)

The reduction of x with respect to the Groebner basis G is -125/48*y^3+41/12*y. We see that any polynomial in x reduces to a polynomial in the monomial basis 1, y, y^2, y^3. So with respect to this particular Groebner basis G every polynomial in the ring has a normal form. The singular.reduce() normalizes polynomials with respect to a Groebner basis. We explain next how this normalization relates to solving a polynomial system.

To set up a multiplication matrix we reduce all multiples of y with respect to this monomial basis.

L = [singular.reduce(y*y^k, G) for k in range(4)]
L

We see that L is [y, y^2, y^3, 8/5*y^2-144/625] The coefficients of these polynomials define a matrix. We first convert the list L to a list of Sage elements.

S = [e.sage() for e in L]
M = [[e.coefficient({y:k}) for k in range(4)] for e in S]
M

The type of the matrix will be over the real numbers.

m = Matrix(RR, 4, 4, M)
m

The matrix m is a companion matrix, its eigenvalues give the values for the y-coordinates of the solutions.

sy = m.eigenvalues()
sy

Observe that the eigenvalues are sorted. To get the corresponding x-coordinates we use the reduction of x with respect to G.

px = singular.reduce(x, G)
spx = px.sage()
sols = [(spx.subs(y=v), v) for v in sy]
sols

We now see the solutions as a list of tuples with values for the x and y coordinates.

Assignments

1. Consider the polynomial system

   \[\begin{split}\left\{ \begin{array}{rcl} x^2 y - y^3 + 4 x + 2 & = & 0 \\ x^3 y + 3 x^2 - 4 y - 9 & = & 0. \end{array} \right.\end{split}\]

   Use a lexicographic Groebner basis to compute the number of solutions.

   What is the number of solutions? Justify your answer.

2. The crank of length \(L\) in a four bar mechanism is centered at \((0, 0)\) and has end points \((L \cos(t), L \sin(t))\) for some angle \(t\). We abbreviate \(\cos(t)\) by \(c\) and \(\sin(t)\) by \(s\). With the introduction of \(c\) and \(s\) we replace the explicit dependency on \(t\) by the implicit equation \(c^2 + s^2 - 1 = 0\). The coordinates \((x,y)\) of the connector point are at the a distance \(r\) of a point with coordinates \((a, 0)\). The connector pont is a distance \(R\) from the end point of the crank. Then we have the system

   \[\begin{split}\left\{ \begin{array}{r} (x-a)^2 + y^2 - r^2 = 0 \\ (x-L c)^2 + (y-L s)^2 - R^2 = 0 \\ c^2 + s^2 - 1 = 0 \end{array} \right.\end{split}\]

   Solve this system for \(x\) and \(y\) symbolically, leaving the six parameters \(L\), \(c\), \(s\), \(a\), \(r\), and \(R\) as symbols

   1. Declare a polynomial ring in eight variables with lexicographic order with rational coefficients. Use the three polynomials to define an ideal in this ring.

   2. Compute a Groebner basis. Examine the structure of the basis. If it does not look triangular to you, then you may have to change the order of the variables as defined in the polynomial ring.

   3. Find explicit formulas for the coordinates \(x\) and \(y\) in terms of the six parameters of the problem. Is there a discriminant expression that is zero for bad values of the parameters?