How to use the package matrixrepresentation in Mathemagix

How to use the package matrixrepresentation in Mathemagix

Luu Ba Thang

GALAAD, INRIA Sophia Antipolis

Email: lbathang@sophia.inria.fr

1.Introduction

Computing the intersection between algebraic varieties is a fundamental task in Computer Aided Geometric Design. Several methods and approaches have been developed for that purpose. Some of them are based on matrix representations of the objects that allow to transform the computation of the intersection locus into generalized eigencomputations . The aim of our works is to compute the matrix representation of these objects and to show that similar algorithms can be implemented even if the matrix representation used are non square matrices. Notice that recent researches, known under the name of the moving lines, moving surfaces, -basis method, have demonstrated that these non square representation matrices are much easier to compute than square representation matrices. In particular, the matrix representation of rational space curves are always non-square matrices

The approach to the curve/curve, curve/surface intersection problem that we develop in our works contain two main parts. The first one is the computation of a matrix representation of the curve and surface from theirs parametrization. This work requires advanced concepts of commutative algebra and algebraic Geometry such as: Elimination theory, free resolution of modules, Koszul complex, approximate complex, Fitting invariant, MacRae invariant, algebraic variety, primary decomposition of ideal,....One of the most effective algorithm for computing them is via Grobner basis of modules, Grobner basis of ideals or via -basis of a set polynomial. All of this computings are symbolic computations.

After combining this matrix representation of the curve and the surface with the parameterization of the curve, the second part consists of a matrix reduction and some eigenvalue computations. In this part, we use the technics of linear algebra such as: Kronecker form of a pencil matrix, Smith form of univariate polynomial matrix, LU-decomposition, QR-decomposition, Sylvester matrix, Bezout matrix,.... As a particularity of our method, these two parts can be performed either by symbolic exact computations or by numerical approached computations. However, it seems that a good combination is to choose a symbolic treatment for the first part, so that the change of representation does not affect the intersection locus, and then using numerical computations, typically LU-decompositions, QR-decomposition and eigenvalues computations, to end the algorithm.

We develop this package in order to solve all problems which we talk above. Hereafter, we introduce some functions for explaining the work which we make.

Mathemagix: http://www.mathemagix.org/www/main/index.en.html

Package: Algebramix (Symbolic Computation).

Packages: Realroot, Numerix, Linalg (Numeric Computation).

2.-basis of a set polynomials

Given the set of homogeneous polynomial of the same degree in with .

By Hilbert-Burch Theorem: Syz() is free and graded -module of rank n. Chosing a basis of .

Definition 1. is called a -basis of .

Remark 2. Given the set of univariate polynomial in with n:=max (deg and . Homogenizing by new variable t which are denoted such that are homogeneous polynomial in K[s,t] of the same degree n. We chose a basis of Syz . Then, the set is called a -basis of the set univariate polynomial.

Function: mubase f gives the -basis of the set of univariate polynomial .

Function: mubase_homogeneous (f,var) gives the -basis of the set of homegenous polynomial with homogeneous variable var.

Example 3. .

Then -basis of is

Mmx]	include "shape/matrixrepresentation/mubase.mmx"

Mmx]	a1:=polynomial(0,0,3);a2:=polynomial(0,3,0,-3); a3:=polynomial(1,1); a4:=polynomial(1,0,1);

Mmx]	f:=[a1,a2,a3,a4]

Mmx]

mubase f

Mmx]	c1:=polynomial(-33,115/2,-49/2,0,1); c2:=polynomial(-36,61,-25,0,1); c3:=polynomial(-8,27/2,-13/2,1);c4:=polynomial(1);

Mmx]	f:=[c1,c2,c3,c4]

Mmx]

mubase f

Mmx]	f0:=QQ['x]<<"x^3+x+1"; f1:=QQ['x]<<"x^4-x+2";f2:=QQ['x]<<"x^6-2*x-1";

Mmx]	f:=[f0,f1,f2]

Mmx]

mubase f

Example 4. .

Then -basis of is:

Mmx]	f0:=QQ['x,'y]<<"x^3+xy^2+y^3"; f1:=QQ['x,'y]<<"x^3-x^2y"; f2:=QQ['x,'y]<<"x^3+x^2y+xy^2"; varf1:=QQ['x,'y]<<"x"; varf2:=QQ['x,'y]<<"y";

Mmx]	f:=[f0,f1,f2]

Mmx]	mubase_homogeneous(f,[varf1,varf2])

3.Matrix representation of parameterized curve and parameterized surface

3.1.Matrix representation of parameterized curve

Let be homogeneous polynomials in of the same degree such that . Consider the regular map

Algebraic curve := Image of which is called a rational curve.

Definition 5. A matrix with entries in is said to be a representation of a parameterized curve if

is generically full rank,

The rank of drops exactly on the curve .

Now, we give some example to compute the matrix representation of a parameterized curve. We focus on curve in and because they are the most interesting cases.

Function: matrix_rep_curve_plane_homogeneous(curplane,varcur,varimplcurve).

Function: matrix_rep_curve_plane(curplane,varimplcurve).

Function: matrix_rep_curve_space_homogeneous(curspace,varcur,varimplcurve).

Function: matrix_rep_curve_space(curspace,varimplcurve).

Example 6. .

Then a matrix representation of is

Mmx]	include "shape/matrixrepresentation/curvesurface.mmx"

Mmx]	R:=QQ['s,'t,'x,'y,'z,'w];

Mmx]	curplane:=[R<<"s^3+t^3",R<<"s^2t",R<<"s^3+st^2+t^3"]

Mmx]	varcur:=[R<<"s",R<<"t"]

Mmx]	varimplcurve:=[R<<"x",R<<"y",R<<"z"]

Mmx]	matrix_rep_curve_plane_homogeneous(curplane,varcur,varimplcurve)

Mmx]	curplane:=[polynomial(1,0,0,1), polynomial(0,0,1),polynomial(1,1,0,1)]

Mmx]	matrix_rep_curve_plane(curplane,varimplcurve)

Example 7. .

Then a matrix representation of is

Mmx]	curspace:=[R<<"3s^2t",R<<"-3s^3+3st^2",R<<"st^2+t^3",R<<"t^3"]

Mmx]	varimplcurspace:=[R<<"x",R<<"y",R<<"z",R<<"w"]

Mmx]	matrix_rep_curve_space_homogeneous(curspace,varcur,varimplcurspace)

Example 8.

Then a matrix representation of is

Mmx]	curspace:=[polynomial(0,0,3),polynomial(0,3,0,-3),polynomial(1,1), polynomial(1)]

Mmx]	matrix_rep_curve_space(curspace,varimplcurspace)

3.2. Matrix representation of parameterized surfaces

Suppose given a parametrization

of a surface such that .

Definition 9. A matrix with entries in is said to be a representation of surface if

is generically full rank,

The rank of drops exactly on the surface .

Example 10. Steine's Roman Surface . Then a matrix representation of is

Function: matrix_rep_surface(surface,varsur,varimplsur).

Mmx]	include "shape/matrixrepresentation/curvesurface.mmx"

Mmx]	R:=QQ['x,'y,'z,'w];

Mmx]	Steine:=[R<<"x^2+y^2+z^2",R<<"xy",R<<"xz",R<<"y*z"]

Mmx]	varimplsur:=[R<<"x",R<<"y",R<<"z",R<<"w"]

Mmx]	varsur:=[R<<"x",R<<"y",R<<"z"]

Mmx]	matrix_rep_surface(Steine,varsur,varimplsur)

Example 11. A sphere :. Then a matrix representation of is

Mmx]	Sphere:=[R<<"x^2+y^2+z^2",R<<"2xz",R<<"2xy",R<<"x^2-y^2-z^2"]

Mmx]	matrix_rep_surface(Sphere,varsur,varimplsur)

Mmx]

4.Polynomial Matrix and Generalized eigenvalues

4.1. Linearization of a polynomial matrix

Let and be two matrices of size . We will call a generalized eigenvalue of A and B a value in the set

Suppose given an -matrix with polynomial entries . It can be equivalently written as a polynomial in with coefficients -matrices with entries in : if then

where .

Definition 12. The generalized companion matrices A, B of the matrix M(t) are the matrices with entries in of size that are given by

where I stands for the identity matrix and stand for the transpose of the matrice .

Function: list_coefficients_matrix M gives the list of matrices: .

Function: firstcompanionmatrix M gives the matrix A.

Function: secondcompanionmatrix M gives the matrix B.

Example 13. M=.

Then and the generalized companion matrices of M are

Mmx]	include "shape/matrixrepresentation/pencilregular.mmx"

Mmx]	M:=matrix(polynomial(0,0,3),polynomial(0,3,0,-3);polynomial(1,1), polynomial(1);polynomial(0,1,2), polynomial(0,1,1))

Mmx]

degree M

Mmx]	list_coefficients_matrix M

Mmx]	firstcompanionmatrix M

Mmx]	secondcompanionmatrix M

Mmx]

4.2.Regular matrix and generalized eigenvalues

Theorem 14.

Definition 15. Pencil matrices is called to be regular pencil matrices if A',B' are square matrices and B' is invertible.

Function: pencilregular(A,B) gives (A',B').

Function: generalized_eigenvalues M gives the list of generalized eigenvalues of M(t).

Example 16. M=.

Regular pencil matrices (A',B') of M and generalized eigenvalues of M are:

Mmx]	include "shape/matrixrepresentation/pencilregular.mmx"

Mmx]	M:=matrix(polynomial(0,0,3),polynomial(0,3,0,-3);polynomial(1,1), polynomial(1); polynomial(0,1,2), polynomial(0,1,1))

Mmx]	A:=firstcompanionmatrix M

Mmx]	B:=secondcompanionmatrix M

Mmx]	pencilregular(A,B)

Mmx]	generalized_eigenvalues M

Mmx]

Example 17. .

Then the generalized eigenvalues of M are:

Mmx]	include "shape/matrixrepresentation/pencilregular.mmx"

Mmx]	M:=matrix(polynomial(0,0,3),polynomial(0,3,0,-3),polynomial(1,2,3,4); polynomial(1,1),polynomial(1), polynomial(-1,3,0,5,6); polynomial(2,1), polynomial(1,3,5), polynomial(0,1,-1,3,4))

Mmx]	generalized_eigenvalues M

Mmx]

5.Parameterized curve/curve intersection

(1)

(2)

: Representation matrix of .

Representation matrix of :

As a consequence of the properties of a representation matrix, we have the following easy property.

Lemma 18. Let , then drops iff the point .

Hereafter, we give some example which are focussed on the intersection of the rational curves in the plane and the intersection of the rational curves in the space .

Function: intersection_curve_plane(curve1,curve2).

Function: intersection_curve_space(curve1,curve2).

Function: intersection_curve_plane_homogeneous(curve1,varcur1,curve2,varcur2).

Function: intersection_curve_space_homogeneous(curve1,varcur1,curve2,varcur2).

Example 19. :

are:

Mmx]	include "shape/matrixrepresentation/curvesurface.mmx"

Mmx]	curve1:=[polynomial(0,0,3),polynomial(0,3,0,-3), polynomial(1,1)]

Mmx]	curve2:=[polynomial(1),polynomial(0,1),polynomial(0,0,1)]

Mmx]	intersection_curve_plane(curve1,curve2)

Mmx]

Figure 1. .

Example 20. .

are:

Mmx]	cur1:=[polynomial(-33,115/2,-49/2,0,1),polynomial(-36,61,-25,0,1), polynomial(-8,27/2,-13/2,1),polynomial(1)]

Mmx]	cur2:=[polynomial(-3,17/2,-11/2,1),polynomial(-6,12,-6,1),polynomial (-38,125/2,-51/2,0,1),polynomial(1)]

Mmx]	intersection_curve_space(cur1,cur2)

Mmx]

Figure 2. Four intersection points : (1,1,0),(0,2,1),(0,3,1) and (-308,-341,-363). The point (-308,-341,-363) doesn't appear in this picture.

Figure 3. The intersection point: (-308,-341,-363)

Example 21.

are:

Mmx]	include "shape/matrixrepresentation/curvesurface.mmx"

Mmx]	cur1:=[QQ['x,'y]<<"x^3+3x^2y+2xy^2+y^3",QQ['x,'y] <<"x^2y+2xy^2+2y^3",QQ['x,'y]<<"xy^2+3y^3", QQ['x,'y]<<"y^3"]

Mmx]	cur2:=[QQ['x,'y]<<"x^4-xy^3+y^4",QQ['x,'y]<<"xy^3+y^4",QQ['x,'y]<<"xy^3+2y^4",QQ['x,'y]<<"y^4"]

Mmx]	varcur:=[QQ['x,'y]<<"x",QQ['x,'y]<<"y"]

Mmx]	intersection_curve_space_homogeneous(cur1,varcur,cur2,varcur)

Mmx]

Figure 4. The intersection points in the affine space: (1,2,3) and (1,1,2).

6.Parameterized curve/surface intersection

matrix representation of the surface .

M(s,t): matrix representation of :

and we have the following easy property:

Lemma 22. Let , then rank of drops iff .

Function: intersection_curve_surface_homogeneous(surface,varsur,curve,varcur)

Example 23. Sphere .

Twisted cubic .

are:

Mmx]	include "shape/matrixrepresentation/curvesurface.mmx"

Mmx]	R:=QQ['s,'t,'x,'y,'z]

Mmx]	sur:=[R<<"x^2+y^2+z^2",R<<"2xz",R<<"2xy",R<<"x^2-y^2-z^2"]

Mmx]	cur:=[R<<"t^3",R<<"st^2",R<<"s^2t",R<<"s^3"]

Mmx]	varsur:=[R<<"x",R<<"y",R<<"z"]

Mmx]	varcur:=[R<<"s",R<<"t"]

Mmx]	intersection_curve_surface_homogeneous(sur,varsur,cur,varcur)

Mmx]

Figure 5. Two intersection points in the affine space : (-0.73735, 0.54368, -0.40089) and (0.73735, 0.54368, 0.40089)

Example 24. Steiner's Roman Surface

Then are:

Mmx]	include "shape/matrixrepresentation/curvesurface.mmx"

Mmx]	steiner:=[R<<"x^2+y^2+z^2", R<<"yz", R<<"xz", R<<"x*y"]

Mmx]	varsteiner:=[R<<"x",R<<"y",R<<"z"]

Mmx]	cur:=[R<<"s^3+st^2+t^3",R<<"s^2t",R<<"s*t^2", R<<"t^3"]

Mmx]	varcur:=[R<<"s", R<<"t"]

Mmx]	intersection_curve_surface_homogeneous(steiner,varsteiner,cur,varcur)

Mmx]

Figure 6. The intersection point in the affine space: (0,0,1), and (0,0,0).

7. Singular points of parameterized plane curve

Let be a rational curve of parameterized by the birational map

Definition 25. A mutiplicity of a point is number of parametervalues (including multiplicities) such that . Denote .

Remark 26. This definition corresponds to the classical definition of multiplicity.

Function: singular_point_plane(curve,varcur,int)

int:=.

Example 27. .

of multiplicity 2 and of multiplicity 3.

Mmx]	include "shape/matrixrepresentation/singular_point.mmx"

Mmx]	R:=QQ['t,'x,'y,'z];

Mmx]	curve1:=[R<<"t^5-2*t^3+2", R<<"t^4+3", R<<"1"]

Mmx]	varcur:=[R<<"t",R<<"x",R<<"y",R<<"z"];

Mmx]	singular_point_plane(curve1,varcur,2)

Mmx]	singular_point_plane(curve1,varcur,3)

Mmx]	singular_point_plane(curve1,varcur,4)

Figure 7. P(2,7) of multiplicity 2 and Q(2,3) of multiplicity 3

Example 28. (Chen and all) .

Mmx]	include "shape/matrixrepresentation/singular_point.mmx"

Mmx]	curve2:=[R<<"t^4-40t^3+40t+1", R<<"t^4+480t^2+1", R<<"t^4+40t^3+480t^2+40t+1"]

Mmx]	singular_point_plane(curve2,varcur,2)

Mmx]	singular_point_plane(curve2,varcur,3)

Example 29. (Busé and D'Andrea).

Mmx]	include "shape/matrixrepresentation/singular_point.mmx"

Mmx]	f:=R<<"t";g:=R<<"t+1";h:=R<<"2t+1";k:=R<<"3t^2+2*t+1";

Mmx]	curve3:=[f^2g^6h^2,f^3h^2k,-g^10]

Mmx]	singular_point_plane(curve3,varcur,2)

Mmx]	singular_point_plane(curve3,varcur,3)

Mmx]	singular_point_plane(curve3,varcur,4)

Mmx]	singular_point_plane(curve3,varcur,6)

Mmx]

8.Conclusion

Introduce new matrix- based representation of rational space curves and parameterized surfaces.
Transfer the solving of the curve/curve, curve/surface intersection problem into the eigenvalues computing problems.
Remove the Kronecker bad blocks of the pencil of matrices in order to extract the regular pencil .
Illustrate the advantages of matrix representation by addressing several impotant problems: the intersection problems and the computation of singularities of rational curves.

9.Some more funtions

9.1.Solve the equation of univariate polynomial

Let be univariate polynomial with the rational coefficients. Funtion solve f is given to find all roots of equation .

Mmx]	include "shape/matrixrepresentation/solvepolynomial.mmx"

Mmx]	f:=QQ['x]<<"x^40+x^3-x^2+x+1"

Mmx]

solve f

Mmx]	g:=polynomial(0,0,1)

Mmx]

solve g^2

9.2.Graded Koszul map

Denote by the set homogemeous polynomial of degree k and then matrix reperesentation of graded Koszul map (in the standard basis)

Function: matrix_multiply_curve(f,var,n,d)

Mmx]	include "shape/matrixrepresentation/curvesurface.mmx"

Mmx]	f:=QQ['s,'t]<<"s^3+2s^2t-t^3"

Mmx]

degree f

Mmx]	var:=[QQ['s,'t]<<"s",QQ['s,'t]<<"t"]

Mmx]	matrix_multiply_curve(f,var,2,3)

Mmx]	matrix_multiply_curve(f,var,3,3)

Denote by the set homogemeous polynomial of degree k and then matrix reperesentation of graded Koszul map (in the standard basis)

Function: matrix_multiply_surface(f,var,n,d)

Mmx]	f:=QQ['s,'t,'u]<<"s^2-2*t^2+u^2"

Mmx]

degree f

Mmx]	var:=[QQ['s,'t,'u]<<"s",QQ['s,'t,'u]<<"t",QQ['s,'t,'u]<<"u"]

Mmx]	matrix_multiply_surface(f,var,2,2)

Mmx]	matrix_multiply_surface(f,var,3,2)

9.3.Extented Euclidean Algorithm

Given two polynomials with . It exits two polynomials of smallest degree such that .

Function: EEA(f,g) give .

Mmx]	include "shape/matrixrepresentation/EEA.mmx"

Mmx]	f:=QQ['x]<<"x^3+3*x-4"

Mmx]	g:=QQ['x]<<"x^5-3x^4+4x-2"

Mmx]

EEA(f,g)

Mmx]

gcd(f,g)

Mmx]	L:=EEA(f,g)

Mmx]	L[1]f+L[2]g

9.4.Another functions

Mmx]	L:=[b1,b2,b3,b4]

Mmx]	delete(L,b1)

Mmx]	delete(L,b4)

Mmx]	delete(L,a)

Mmx]	L:=[b1,b1,b2,b2,b3,b4,b4,b4]

Mmx]	delete_dif L

Mmx]

10.Some future works

Using subresultant of -basis to solve some problems of rational space curves.
Intersection multiplicity problem between two rational curves.
Intersection parameterized surface/surface.
Singular points of parameterized surfaces.
Develop the implementation of the software Mathemagix to solve the above purposes.

11. References

1. N. Song, R. Goldman, -bases for polynomial systems in one variable, Comput. Aided Geom. Design 26 (2) (2009) 217–230.

2. T. Luu Ba, L. Busé, B. Mourrain, Curve/surface intersection problem by means of matrix representations, in: SNC '09: Proceedings of the 2009 conference on Symbolic numeric computation, ACM, New York, NY, USA, 2009, pp. 71–78.

3. L. Busé and M. Chardin, Implicitizing rational hypersufaces using approximation complexes. J. Symbolic Comput. 40(4-5): 1150-1168, 2005.

4. L. Busé and T. Luu Ba, Matrix-based implicit representations of rational algebraic curves and applications. To appear Computed Aided Geometric Design 2010.

5. F. Chen, W. Wang and Y. Liu, Computing singular points of plane rational curves. J. Symbolic Comput. 43 , 92-117, 2008.

6. The software Axel (http://axel.inria.fr) and Mathemagix (http://www.mathemagix.org).