Computational Geometry and Computer Algebra Sylvain Lazard INRIA - - PowerPoint PPT Presentation

Computational Geometry Computer Algebra

and

• C. Borcea, Rider Univ., New Jersey
• H. Br¨
• nnimann, Poly. Univ., New York
• O. Devillers, Inria Sophia-Antipolis
• L. Dupont, H. Everett, X. Goaoc, S. Petitjean, Loria - Inria Nancy
• D. Lazard, M. Safey El Din, UPMC Paris
• F. Sottile, Texas A&M Univ., Texas

Sylvain Lazard

INRIA Nancy Joint work with

Context

Development of efficient algorithms and data structures for solving geometric problems

Computational Geometry

(since the 1970’s)

Historically: Problems with points, segments, polygons, meshes in 2, 3 and d dimensions

Context

Development of efficient algorithms and data structures for solving geometric problems

Closest pair of points Convex hull Delaunay triangulation Voronoi diagram

Computational Geometry

(since the 1970’s)

Historically: Problems with points, segments, polygons, meshes in 2, 3 and d dimensions

Context

Development of efficient algorithms and data structures for solving geometric problems

Closest pair of points Convex hull Delaunay triangulation Voronoi diagram

Computational Geometry

(since the 1970’s)

Historically: Problems with points, segments, polygons, meshes in 2, 3 and d dimensions

Context

Development of efficient algorithms and data structures for solving geometric problems

Closest pair of points Convex hull Delaunay triangulation Voronoi diagram

Computational Geometry

(since the 1970’s)

Historically: Problems with points, segments, polygons, meshes in 2, 3 and d dimensions

Context

Development of efficient algorithms and data structures for solving geometric problems

Closest pair of points Convex hull Delaunay triangulation Voronoi diagram

Computational Geometry

(since the 1970’s)

Historically: Problems with points, segments, polygons, meshes in 2, 3 and d dimensions

Context

Development of efficient algorithms and data structures for solving geometric problems

Closest pair of points Convex hull Delaunay triangulation Voronoi diagram

Computational Geometry

(since the 1970’s)

Historically: Problems with points, segments, polygons, meshes in 2, 3 and d dimensions

Context

Development of efficient algorithms and data structures for solving geometric problems

Computational Geometry

(since the 1970’s)

Historically: Problems with points, segments, polygons, meshes in 2, 3 and d dimensions

Line segment intersection Linear programming Closest pair of points Convex hull Delaunay triangulation Voronoi diagram

Context

Development of efficient algorithms and data structures for solving geometric problems

Computational Geometry

(since the 1970’s)

Historically: Problems with points, segments, polygons, meshes in 2, 3 and d dimensions

Line segment intersection Linear programming Closest pair of points Convex hull Delaunay triangulation Voronoi diagram

Context

Development of efficient algorithms and data structures for solving geometric problems

Computational Geometry

(since the 1970’s)

Historically: Problems with points, segments, polygons, meshes in 2, 3 and d dimensions

Line segment intersection Linear programming Closest pair of points Convex hull Delaunay triangulation Voronoi diagram

Context

Development of efficient algorithms and data structures for solving geometric problems

Computational Geometry

(since the 1970’s)

Historically: Problems with points, segments, polygons, meshes in 2, 3 and d dimensions

Line segment intersection Linear programming Closest pair of points Convex hull Delaunay triangulation Voronoi diagram Meshing

Context

Development of efficient algorithms and data structures for solving geometric problems

More recently: curves and surfaces in 2D and 3D Non-Linear Computational Geometry

❀

Connexions with Computer algebra, Symbolic comput. Computational Geometry

(since the 1970’s)

Historically: Problems with points, segments, polygons, meshes in 2, 3 and d dimensions

Context

Manipulation of curved objects

CAD/CAGD, solid modeling, computer graphics

Context

Manipulation of curved objects

CAD/CAGD, solid modeling, computer graphics Most software (even commercial) is not robust

Context

Manipulation of curved objects

CAD/CAGD, solid modeling, computer graphics Most software (even commercial) is not robust They choke, crash, loop or return wrong answers

• n (near)-degenerate instances

(aligned objects, tangent primitives, multiple intersections)

Robustness and exact algorithms

Geometric decisions: Is R to the left, to the right or on line PQ? P Q R

Robustness and exact algorithms

Geometric decisions: Is R to the left, to the right or on line PQ? P Q R Fixed-precision floating-point arithmetic ❀ Errors in the geometric decisions ❀ Major errors & crashes

Robustness and exact algorithms

Geometric decisions: Is R to the left, to the right or on line PQ? P Q R Fixed-precision floating-point arithmetic ❀ Errors in the geometric decisions

Ex: P = (0.5000029, 0.5000027), Q = (12, 12), R = (24, 24)

(xQ − xP )(yR − yP )− (xR − xP )(yQ − yP ) = { Float: +3.05 . . . 10−5 ❀ R is left of PQ Exact: −2.410−6 ❀ R is right of PQ

❀ Major errors & crashes

Robustness and exact algorithms

Geometric decisions: Is R to the left, to the right or on line PQ? P Q R Fixed-precision floating-point arithmetic ❀ Errors in the geometric decisions Ex: Convex hull of points ❀ Major errors & crashes Point missed Not convex

[Interstices]

Robustness and exact algorithms

Geometric decisions: Is R to the left, to the right or on line PQ? P Q R Fixed-precision floating-point arithmetic ❀ Errors in the geometric decisions Ex: Convex hull of points ❀ Major errors & crashes Point missed Not convex

[Interstices]

Robustness and exact algorithms

Paradigm of exact geometric computing Implement geometric decisions exactly Treat all degenerate situations (Geometric constructions may be approximate)

Robustness and exact algorithms

Implement with exact arithmetic Paradigm of exact geometric computing Implement geometric decisions exactly Treat all degenerate situations Solve geometric problems exactly (Geometric constructions may be approximate) In other words use speed-up for efficiency, e.g. interval arithmetic

Robustness and exact algorithms

Degeneracies may be hard to characterize/handle Exact computation may be too slow General philosophy for the design of geometric algorithms Avoid (as much as possible) algebraic computations and the appearance of algebraic numbers due to high-degree algebraic numbers Bad news Rely (as much as possible) on the geometry

Computer algebra in computational geometry

Off line computation/analysis

Study of base cases Characterization of degeneracies

Rely as much as possible on the geometry

❀

Computer algebra & theorems in geometry Algebraic computations in geometric algorithms

Topology of algebraic curves Algebraic kernel in comp. geom. library CGAL Intersection of quadrics

Avoid as much as possible algebraic computations but still some computation is needed. . . Theorems

Plan

I.a. Properties of tangent lines in 3D

• II. Avoiding algebraic numbers

Intersections of two quadrics

• I. Computer algebra & theorems in geometry

I.b. Voronoi diagrams of three lines

Properties of tangent lines in 3D

Key structure for 3D visibility: Visibility complex [PV 96, DDP 02] Partition into connected components of segments that touch the same objects set of maximal free line segments

light source penumbra umbra

Properties of tangent lines in 3D

Key structure for 3D visibility: Visibility complex [PV 96, DDP 02] Partition into connected components of segments that touch the same objects Vertices: set of maximal free line segments

• max. free line segment tangent to 4 objects

light source penumbra umbra

Properties of tangent lines in 3D

Problem 1. What is the maximum number of lines tangent to 4 triangles?

Properties of tangent lines in 3D

Problem 1. What is the maximum number of lines tangent to 4 triangles? 10? 300?

Properties of tangent lines in 3D

Problem 1. What is the maximum number of lines tangent to 4 triangles?

4 “fat” triangles with 40 tangents

Properties of tangent lines in 3D

Problem 1. What is the maximum number of lines tangent to 4 triangles?

4 “fat” triangles with 40 tangents

Lower bound: 62 Upper bound: 162 (naive: 4 · 34 = 324)

Techniques: geometry of quadrics, combinatorics

[Br¨

• nnimann, Devillers, Lazard, Sottile, Discrete & Comput. Geom., 2007]

experiments,

Properties of tangent lines in 3D

Problem 1. What is the maximum number of lines tangent to 4 triangles?

Constant size problem ❀

Partition the configuration space (R30) into components that have the same number of tangents Compute the number of tangents for one instance in each component Can be solved in theory, but in practice. . .

(Degeneracies may be hard to characterize) Do 4 balls admit an infinite number of tangents? Problem 2.

Properties of tangent lines in 3D

(Degeneracies may be hard to characterize) Do 4 balls admit an infinite number of tangents? Test if an algebraic system of degree 12 has infinitely many real solutions

[MPT01]

x4 + 17xy3 + . . . = 0 5y3 + 3x2y + . . . = 0

{

Problem 2.

Properties of tangent lines in 3D

(Degeneracies may be hard to characterize) Do 4 balls admit an infinite number of tangents? Test if an algebraic system of degree 12 has infinitely many real solutions

[MPT01]

x4 + 17xy3 + . . . = 0 5y3 + 3x2y + . . . = 0

{

Problem 2. Characterize the configurations of 4 balls that admit infinitely many tangents

Properties of tangent lines in 3D

Lines tangent to balls

Number of tangents to 4 balls: 12 or ∞

[MPT01, DMPT03] [ST02]

Configurations with ∞ many tangents: 4 unit balls: aligned centers

[MPT01]

4 arbitrary balls?

[Borcea, Goaoc, Lazard, Petitjean, Discrete & Comput. Geom., 2006]

4 balls admit infinitely many tangents iff aligned centers and at least 1 common tangent

Lines tangent to balls

Sketch of proof Line represented by (p, v) with p · v = 0 O p − → v

Lines tangent to balls

(1/2)

Sketch of proof Line represented by (p, v) with p · v = 0 Tangency condition with spheres centered at ci of radii ri p2 = r2

1

(with c1 = O) O p − → v

[MPT01]

Lines tangent to balls

2|v|2 @ cT

2

cT

3

cT

4

1 A p = |v|2 @ |c2|2 + r2

1 − r2 2

|c3|2 + r2

1 − r2 3

|c4|2 + r2

1 − r2 4

1 A − @ (c2 · v)2 (c3 · v)2 (c4 · v)2 1 A Φ2(v)

{

Φ0

{

M

{

(1/2)

Sketch of proof Line represented by (p, v) with p · v = 0 Tangency condition with spheres centered at ci of radii ri p2 = r2

1

(with c1 = O) O p − → v

[MPT01]

Lines tangent to balls

2|v|2 @ cT

2

cT

3

cT

4

1 A p = |v|2 @ |c2|2 + r2

1 − r2 2

|c3|2 + r2

1 − r2 3

|c4|2 + r2

1 − r2 4

1 A − @ (c2 · v)2 (c3 · v)2 (c4 · v)2 1 A Φ2(v)

{

Φ0

{

M

{

(1/2)

Idea: For each variable, eliminate all other variables Infinite number of solutions ⇔ all the coefficients Deduce the geometric characterization

• f the uivariate polynomials are 0

Sketch of proof Line represented by (p, v) with p · v = 0 Tangency condition with spheres centered at ci of radii ri p2 = r2

1

(with c1 = O) O p − → v

[MPT01]

Lines tangent to balls

2|v|2 @ cT

2

cT

3

cT

4

1 A p = |v|2 @ |c2|2 + r2

1 − r2 2

|c3|2 + r2

1 − r2 3

|c4|2 + r2

1 − r2 4

1 A − @ (c2 · v)2 (c3 · v)2 (c4 · v)2 1 A Φ2(v)

{

Φ0

{

M

{

(1/2)

Idea: For each variable, eliminate all other variables Infinite number of solutions ⇔ all the coefficients Deduce the geometric characterization In practice, we did not succeed. . .

• f the uivariate polynomials are 0

Sketch of proof Line represented by (p, v) with p · v = 0 Tangency condition with spheres centered at ci of radii ri p2 = r2

1

(with c1 = O) O p − → v Proof, case 1: Suppose non-coplanar centers (M inversible) Prove that the 4 balls admit finitely many tangents

[MPT01]

Lines tangent to balls

2|v|2 @ cT

2

cT

3

cT

4

1 A p = |v|2 @ |c2|2 + r2

1 − r2 2

|c3|2 + r2

1 − r2 3

|c4|2 + r2

1 − r2 4

1 A − @ (c2 · v)2 (c3 · v)2 (c4 · v)2 1 A Φ2(v)

{

Φ0

{

M

{

(1/2)

Sketch of proof Line represented by (p, v) with p · v = 0 Tangency condition with spheres centered at ci of radii ri p2 = r2

1

(with c1 = O) O p − → v Proof, case 1: Suppose non-coplanar centers (M inversible) Express p in terms of v: Substitute p in p · v = 0 and p2 = r2

1

p = M −1

1 2Φ0 + 1 2|v|2 Φ2(v)

• Prove that the 4 balls admit finitely many tangents

[MPT01]

Lines tangent to balls

2|v|2 @ cT

2

cT

3

cT

4

1 A p = |v|2 @ |c2|2 + r2

1 − r2 2

|c3|2 + r2

1 − r2 3

|c4|2 + r2

1 − r2 4

1 A − @ (c2 · v)2 (c3 · v)2 (c4 · v)2 1 A Φ2(v)

{

Φ0

{

M

{

(1/2)

• M −1

Φ2(v) + |v|2Φ0

• · v = 0
• M −1

Φ2(v) + |v|2Φ0

• 2 = 4r2

1|v|4

(cubic) (quartic)

{

Lines tangent to balls

Sketch of proof O p − → v

(2/2)

• M −1

Φ2(v) + |v|2Φ0

• · v = 0
• M −1

Φ2(v) + |v|2Φ0

• 2 = 4r2

1|v|4

(cubic) (quartic)

Suppose for a contradiction that there are ∞ many solutions

The solutions form a curve Γ in P2(R)

{

Lines tangent to balls

Sketch of proof O p − → v

(2/2)

v ∈ P2(R)

• M −1

Φ2(v) + |v|2Φ0

• · v = 0
• M −1

Φ2(v) + |v|2Φ0

• 2 = 4r2

1|v|4

(cubic) (quartic)

Suppose for a contradiction that there are ∞ many solutions

The solutions form a curve Γ in P2(R) In P2(C), Γ intersect any curve and, in particular,

{

Lines tangent to balls

Sketch of proof O p − → v

(2/2)

v ∈ P2(R) the conic |v|2 = 0

• M −1

Φ2(v) + |v|2Φ0

• · v = 0
• M −1

Φ2(v) + |v|2Φ0

• 2 = 4r2

1|v|4

(cubic) (quartic)

Suppose for a contradiction that there are ∞ many solutions

The solutions form a curve Γ in P2(R) In P2(C), Γ intersect any curve and, in particular,

{

M −1Φ2(v) ≡ v ❀ v is real and |v|2 = 0 ❀ v = 0, contradiction

Lines tangent to balls

Sketch of proof O p − → v

(2/2)

v ∈ P2(R) the conic |v|2 = 0

M −1Φ2(v) is on the conic {|x|2 = 0}

• M −1Φ2(v)
• · v = 0
• M −1Φ2(v)
• 2 = 0

{

M −1Φ2(v) is on the tangent to this conic at v

Properties of tangent lines in 3D

Problem 3. Conjecture: 4 disjoint unit balls (nonaligned centers) admit at most 8 tangents

Properties of tangent lines in 3D

Problem 3. Conjecture: 4 disjoint unit balls (nonaligned centers) Constant size problem admit at most 8 tangents 5 variables: the line 6 parameters: the ball centers 7 inequalities: the balls are disjoint & not aligned 5 equations: tangency and line conditions

Properties of tangent lines in 3D

Problem 3. Conjecture: 4 disjoint unit balls (nonaligned centers) Constant size problem admit at most 8 tangents 5 variables: the line 6 parameters: the ball centers 7 inequalities: the balls are disjoint & not aligned Tried a few years ago and failed, but . . . 5 equations: tangency and line conditions

Properties of tangent lines in 3D

Problem 3. Conjecture: 4 disjoint unit balls (nonaligned centers) Constant size problem admit at most 8 tangents 5 variables: the line 6 parameters: the ball centers 7 inequalities: the balls are disjoint & not aligned Tried a few years ago and failed, but . . . 2

` M −1 ` Φ2(v) + |v|2Φ0 ´´ · v = 0 ˛ ˛M −1 ` Φ2(v) + |v|2Φ0 ´˛ ˛2 = 4|v|4

{

5 equations: tangency and line conditions 2

Properties of tangent lines in 3D

Problem 3. Conjecture: 4 disjoint unit balls (nonaligned centers) Constant size problem admit at most 8 tangents 5 variables: the line 6 parameters: the ball centers 7 inequalities: the balls are disjoint & not aligned Tried a few years ago and failed, but . . . 2

` M −1 ` Φ2(v) + |v|2Φ0 ´´ · v = 0 ˛ ˛M −1 ` Φ2(v) + |v|2Φ0 ´˛ ˛2 = 4|v|4

{

5 equations: tangency and line conditions 2 degree 3/7, 68 monomials degree 4/12, 763 monomials

Properties of tangent lines in 3D

Problem 3. Conjecture: 4 disjoint unit balls (nonaligned centers) Constant size problem admit at most 8 tangents 5 variables: the line 6 parameters: the ball centers 7 inequalities: the balls are disjoint & not aligned Tried a few years ago and failed, but . . . 2

` M −1 ` Φ2(v) + |v|2Φ0 ´´ · v = 0 ˛ ˛M −1 ` Φ2(v) + |v|2Φ0 ´˛ ˛2 = 4|v|4

{

5 equations: tangency and line conditions 2 degree 3/7, 68 monomials degree 4/12, 763 monomials Discriminant variety & a point per connected component of the complement

Properties of tangent lines in 3D

Problem 3. Conjecture: 4 disjoint unit balls (nonaligned centers) Constant size problem admit at most 8 tangents 5 variables: the line 6 parameters: the ball centers 7 inequalities: the balls are disjoint & not aligned Tried a few years ago and failed, but . . . 2

` M −1 ` Φ2(v) + |v|2Φ0 ´´ · v = 0 ˛ ˛M −1 ` Φ2(v) + |v|2Φ0 ´˛ ˛2 = 4|v|4

{

5 equations: tangency and line conditions 2 degree 3/7, 68 monomials degree 4/12, 763 monomials Discriminant variety & a point per connected component of the complement Still too big. . . so far

Plan

I.a. Properties of tangent lines in 3D

• II. Avoiding algebraic numbers

Intersections of two quadrics

• I. Computer algebra & theorems in geometry

I.b. Voronoi diagrams of three lines

Voronoi diagrams

Voronoi diagrams, Dirichlet tessellations, Medial axes In 2D: Well understood for points and polygons In 3D: Well understood for points

Voronoi diagrams

Voronoi diagrams, Dirichlet tessellations, Medial axes In 2D: Well understood for points and polygons In 3D: Well understood for points What about Voronoi diagrams of polyhedra?

• f lines?

Voronoi diagrams of lines

Complexity: n lines: Ω(n2), O(n3+ǫ)

[Sharir94]

n lines with c possible directions: O(n2+ǫ) (conjectured near quadratic)

[KS03]

Algorithms: Exact: naive

(compute, for each line, a cell of an arrangement of quadrics)

Approximation: many Applications: Medial axis of polyhedra

Voronoi diagrams of 3 lines

The topology of the Voronoi diagram is invariant Trisector: 4 infinite branches of a smooth quartic 2D cells: 2 connected components of a hyperb. paraboloid bounded by 3 branches and 1 branch of the trisector, resp. 3 lines in general position

(not parallel to a common plane and pairwise skew)

• r a cubic and a line

[Everett, Lazard, Lazard, Safey El Din, ACM Symp. on Comp. Geom., 2007]

Voronoi diagrams of 3 lines

• Th. All 4 branches are monotonic in Y
• Th. Localizing and ordering points on a branch

can be done by evaluating the sign of linear forms with rational coefficients

[Everett, Lazard, Lazard, Safey El Din, ACM Symp. on Comp. Geom., 2007]

Complete characterization of the Voronoi of three lines ❀ fundamental properties & point-location algorithms

Voronoi diagrams of 3 lines

Idea of the proof

• Trisector of 3 lines in general position:

non-singular with 4 simple real points at ∞ ❀ homeomorphic to 4 lines with no intersection point

• Study of an example
• Topology of the Voronoi diagram is

invariant by continuous deformation

• ver the set of all triplets of 3 lines in gen. position

which is connected

Voronoi diagrams of 3 lines

Idea of the proof

• Trisector of 3 lines in general position:

non-singular ❀ homeomorphic to 4 lines with no intersection point

• Study of an example
• Topology of the Voronoi diagram is

invariant by continuous deformation

• ver the set of all triplets of 3 lines in gen. position

which is connected with 4 simple real points at ∞

Voronoi diagrams of 3 lines

Bisector of two skew lines: Hyperbolic paraboloid (HP)

Voronoi diagrams of 3 lines

Trisector of three lines Intersection of 2 HP bisectors: H1,2 ∩ H1,3 Intersection of 2 quadrics: smooth quartic, nodal quartic, a cubic and a line, 2 conics, 4 lines, etc. Idea of the proof Determine the type of H1,2 ∩ H1,3 Deduce that the trisector is non-singular

Voronoi diagrams of 3 lines

Type of the intersection of two quadrics H1,2, H1,3

Voronoi diagrams of 3 lines

Type of the intersection of two quadrics H1,2, H1,3 Q1,2, Q1,3: 4 × 4 matrices associated with H1,2, H1,3 Pencil of quadrics/matrices: {P(λ) = Q1,2 + λQ1,3 | λ ∈ R}

Voronoi diagrams of 3 lines

Type of the intersection of two quadrics H1,2, H1,3 Q1,2, Q1,3: 4 × 4 matrices associated with H1,2, H1,3 Pencil of quadrics/matrices: {P(λ) = Q1,2 + λQ1,3 | λ ∈ R} Characteristic polynomial D(λ)= det P(λ), degree 4 in λ

Voronoi diagrams of 3 lines

Type of the intersection of two quadrics H1,2, H1,3 Q1,2, Q1,3: 4 × 4 matrices associated with H1,2, H1,3 Pencil of quadrics/matrices: {P(λ) = Q1,2 + λQ1,3 | λ ∈ R} Characteristic polynomial D(λ)= det P(λ), degree 4 in λ

[Segre 1883]

H1,2 ∩ H1,3 is a smooth quartic or ∅ When D(λ) has no multiple root:

Voronoi diagrams of 3 lines

Type of the intersection of two quadrics H1,2, H1,3 Q1,2, Q1,3: 4 × 4 matrices associated with H1,2, H1,3 Pencil of quadrics/matrices: {P(λ) = Q1,2 + λQ1,3 | λ ∈ R} Characteristic polynomial D(λ)= det P(λ), degree 4 in λ

[Segre 1883]

The real type of H1,2 ∩ H1,3 follows from: number of multiples roots, λi, of D(λ), number of positive/negative eigenvalues of P(λi)

[DLLP08]

H1,2 ∩ H1,3 is a smooth quartic or ∅ When D(λ) has no multiple root: When D(λ) has multiple root(s):

Voronoi diagrams of 3 lines

Proof. We compute the characteristic polynomial D(λ) discriminant ∆ its We determine the type of H1,2 ∩ H1,3 when ∆ = 0 (when ∆ = 0, H1,2 ∩ H1,3 is a non-singular quartic)

Voronoi diagrams of 3 lines

The setting Any two skew lines can be moved into the lines Z = 1, Y = a X and Z = −1, Y = −a X The 3 lines are not parallel to a common plane ⇒ Third line: defined by point (x, y, 0) and vector (α, β, 1)

Y X O Z (x, y, 0) (α, β, 1) ℓ3 ℓ2 ℓ1 Y = aX Z = 1 Y = −ax Z = −1

{ {

D(λ) = ` α2 + β2 + 1 ´ a2λ4 − 2 a ` 2 aβ2 + ayβ + aα x − β α + 2 a + 2 aα2 − β α a2´ λ3 +(β2+6 a2β2−2 β xa3−6 β α a3+6 yβ a2−6 aβ α−2 aβ x+6 α xa2 + y2a2 − 2 aα y + x2a2 − 2 yα a3 + 6 a2α2 + a4α2 + 4 a2)λ2 − 2 ` xa − ya2 − 2 β a2 − β + 2 aα + α a3´ (xa − y − β + aα) λ + ` 1 + a2´ (xa − y − β + aα)2

Voronoi diagrams of 3 lines

∆ =16 a4 (a x − y − β + a α)2 (y + a x − a α − β)2 ∗ gros facteur

By our general position assumption, 16 a4 (a x − y − β + a α)2 (y + a x − a α − β)2 > 0

Voronoi diagrams of 3 lines

Key problem: determine the type of the trisector when gros facteur = 8 a8α4y2 + 7 a4β2x4 − 4 aβ3x + 16 a8β4x4 + 32 a4α2y2 +

2 a6α2β4x2 + 38 a8α2x2 + 2 y4β2a4α2 + 44 a8α2β2x2 + 2 a6α4y2β2 + a2α2β4y2 − 8 y4a4 + 3 a8α4β2 − 8 a8α4 + 20 y2β2a2 + β6x2a2 + 20 a4β2α2 + 38 y2β2a4 + 2 α2y4a6x2 + 8 a4β2 + 32 a4y2 − 8 a4x2 − 8 a4β4 + 28 α2y4a5xβ + 10 a7α3β2y+28 a7β2x4yα−4 aβ5x+48 β xa5 +16 β xa3 +32 yα a5 −8 β5xa3 + 48 β3xa5 +32 a7β3x+8 β3xa3 +32 β a7x+a12α6 −28 β α3a8yx+20 β2α2a8 + 8 β2α a7y − 4 β2α a5y + 14 β4α a3y + 10 a5α2β3x + 2 a7α2β3x + 32 a5α3y + 2 a5α3β2y−4 a7α2β x+8 a5α2β x−4 β3α3a6xy+48 a7α3y+x4a4+10 aβ4α y+ 8 aβ2α y + 28 a5β4x2α y − 8 a8x4 + 48 ya7α − 4 a6β x3α y3 − 2 a6β x5α y + 2 a4β4x2y2 + 22 a6α2y2β2 + 32 a2α2y2β2 − 8 a10α2 + 20 a6α2x2 + 16 a4α2 + 20 a6y2β2 + 38 a6β2α2 + 16 a6α2y2 + 16 a6β2x2 + 14 y4a5xβ + 16 a10x2 − 8 a6y4 +8 a6α4y2 +22 a6α2β2x2 +44 a4α2y2β2 +3 a4α2β4 +16 a4 +32 x2a8 + 8 a6β4x2 + 8 a8α2 + a8y4 − 4 a5x5β + β6 − 28 a4β3yα x + 8 a5α3y3x2 + 3 β4y2a2 + 8 β4a4x2 − 12 a6β3xα y + 24 β a6yα x + 32 a6β2 + 8 a6y2 + 8 a6x2 + β4 −8 β4a2 +6 β2x2a2 +16 α2y2a2 −8 y5a5α+3 y4a6x2 +8 a6β2x4 +32 a6α2 − 40 a7β2x2yα + 32 β2x2a8 − 8 a2β2 + 32 a6 − 8 a5β5x3 + 16 y2a2

Voronoi diagrams of 3 lines

−16 β xα ya2 − 8 a8y2 + 7 β4x2a2 + 16 a8β2 + 16 a3β4x2α y + a12α4 + 28 a5β2x2α y−56 a6β3x3α y+54 a4α2y2x2β2 +16 a8 −8 a9α3y3 −8 α3y5a5 − 40 a5α2y2β x − 48 a4α3y3xβ + 44 a6y2α2x2 − 12 y3xa4β α − 4 yα a7x2 + 8 a3yx2α + 8 a3y2xβ − 4 y2xa5β + 32 a3α y3 + 48 y3α a5 − 4 x3a3β + 20 x2a4y2 + 8 x3a5β + 10 a9α3yx2 + 14 a9α4β x + 10 y3β2a5α − 4 a3β5x3 − 8 a10α4 + 2 β3y2a5x + 28 a7α2y2xβ + 8 a7y3α + 8 yα3a9 + 48 a7β3x3 + 8 a6α2y4 +44 a6β2x2y2 +22 a8α2y2x2 −8 a9α5y +8 yx2a9α−32 x2a9yβ2α+ 32 a10β2α2x2 + 20 a10α2x2 + 3 a10α4x2 + 20 y2x2a8 + 32 x2a8y2β2 + 3 y4β2a4 +8 a8β2x4 −32 a5β2x2y3α+2 a5x2y3α+12 a5β3x3y2 +32 a9β3x3 + 32 a9x3β+2 a9α2x3β+54 a8β2x2α2y2−48 a8β3x3α y+6 a8α4y4+7 a8α2y4− 28 a8α3y3β x + 28 a9α2y2β x + β2a10α4x2 + 2 β2a8α2x4 − 2 a10α5yxβ − 24 a10α3yxβ+10 a11α4xβ−4 a9y3α−24 a8β xα y3+8 a9y2β x−12 a8x3β α y+ 16 a9α4xβ y2 + 8 a9x3α2β3 − 32 a7x2α3β2y − 20 a9β2x2α3y − 4 a9α5y3 + 16 a9yα + 12 a7x3α2β3 + 28 a7α4xβ y2 − 24 a4β x3α y + 32 a4y2α2x2 + 2 a8x2α4y2 + 8 α2xa11β + 80 a6x2α2y2β2 − 8 a3β3x3 − 4 α3x3a8β y + 32 a3α3y3 + 16 β2x2a10 + 10 a5x4yα − 20 a3α2y2xβ3 − 12 a6α3yxβ+

Voronoi diagrams of 3 lines

8 a3α3y3β2 + 10 y2β3a3x + 2 y3β2a3α − 16 β xa10yα + 10 β xa7y4 + 12 y3a5α3β2 − 8 y3a7α5 − 4 y5a7α3 + 6 y2α2a10 − 4 a11α5y + a10α6y2 + 7 y2α4a10 − 32 β3xa5α2y2 − 28 y3xa6β α − 56 y3xa6β α3 − 20 a7β2x2α y3 + 8 a7β3x3y2 + 16 a7α2y4xβ − 8 a6x4 − 4 a7y5α + 2 a7x3y2β + 10 a7x2α y3 + 12 a7x2α3y3 − 32 a7x3α2y2β − 4 yα3a11 + 16 y4a6α4 − 20 a5α2y2x3β + 3 a6y2x4+3 a8x4α2−8 a7x5β−8 a7β3x5+6 a4β4x4+48 a5α3y3+16 a4α4y4+ 8 a4y4α2 + 38 a6y2x2 + 16 yα a5β2x4 + 48 β x3a7 + 16 a6β4x4 + a6β2x6 − 2 β5α a2xy +y4β2a6x2 +10 a5x3y2β +14 a7x4yα+2 a6β2x4y2 +a6α2y2x4 − 4 a4β3xα y3 −28 a4β3x3α y −28 a6β x3α y −4 β3x5a5 +22 a4y2β2x2 +y6a6 + α2y6a6 − 2 β xα y5a6 + x6a6 + 10 β x3a7α2 + 2 yα3a7x2 − 32 a3α2y2β x + 28 a3β2x2α y − 24 a2β3yα x

Degree 18 in 5 variables with 253 monomials

is zero

Voronoi diagrams of 3 lines

Technique: Prove that the gros facteur is never negative Manipulation of Gr¨

• bner bases

❀ gros facteur = 0 only if (y + aα)(ax + β) = 0 ❀ gros facteur = 0 iff the trisector is a cubic and line heavy use of computer algebra tools Classification of quadric intersections

(RAGLib) (GB)

• without intersection

Voronoi diagrams of 3 lines

• Lemma. The gros facteur (and ∆) are never negative

Voronoi diagrams of 3 lines

• Lemma. The gros facteur (and ∆) are never negative
• Proof. Compute at least a point in every connected

component of the semi-algebraic set ∆(a, x, y, α, β) < 0

Voronoi diagrams of 3 lines

• Lemma. The gros facteur (and ∆) are never negative
• Proof. Compute at least a point in every connected

component of the semi-algebraic set ∆(a, x, y, α, β) < 0 Major problem in computer algebra Recent advances provide practical solutions

Voronoi diagrams of 3 lines

• Lemma. The gros facteur (and ∆) are never negative
• Proof. Compute at least a point in every connected

component of the semi-algebraic set ∆(a, x, y, α, β) < 0 Major problem in computer algebra RAGLib

[Safey El Din] based on FGb/RS [Faug` ere, Rouillier]

On {(a, x, y, α, β) | ∆ < 0}: returns ∅

(in 10 hours)

Recent advances provide practical solutions

Voronoi diagrams of 3 lines

• Lemma. The gros facteur (and ∆) are never negative
• Proof. Compute at least a point in every connected

component of the semi-algebraic set ∆(a, x, y, α, β) < 0 Major problem in computer algebra RAGLib

[Safey El Din] based on FGb/RS [Faug` ere, Rouillier]

On {(a, x, y, α, β) | ∆ < 0}: returns ∅

(in 10 hours)

∆ = 0 ⇔ the gros facteur and all its partial derivatives are 0 Recent advances provide practical solutions

RAGLib

Finding (at least) a point in every connected component

• f a semi-algebraic set ∆(a, x, y, α, β) < 0

z = ∆(a, x, y, α, β)

RAGLib

Finding (at least) a point in every connected component

• f a semi-algebraic set ∆(a, x, y, α, β) < 0
• Compute all the generalized critical values of ∆
• Compute the largest negative critical value, zc

z = ∆(a, x, y, α, β) zc

RAGLib

Finding (at least) a point in every connected component

• f a semi-algebraic set ∆(a, x, y, α, β) < 0
• Compute all the generalized critical values of ∆
• Compute the largest negative critical value, zc
• Pick a value ˜

z in (zc, 0)

z = ∆(a, x, y, α, β) ˜ z zc

RAGLib

Finding (at least) a point in every connected component

• f a semi-algebraic set ∆(a, x, y, α, β) < 0
• Compute all the generalized critical values of ∆
• Compute the largest negative critical value, zc
• Pick a value ˜

z in (zc, 0)

• Compute (at least) a point

in every connected component of the algebraic set ∆(a, x, y, α, β) = ˜ z

z = ∆(a, x, y, α, β) ˜ z zc

Finding (at least) a point in every connected component

RAGLib

• f the algebraic set ∆(a, x, y, α, β) = ˜

z Solve X × grad(∆)(X) = 0

with X = (a, x, y, α, β)

❀ At least a point in every cc of the semi-algebraic set ∆ < 0 ∆(X) = ˜ z

{

O

Voronoi diagrams of 3 lines

Main Lemma ∆ = 0 only if (y + a α) (a x + β) = 0

Recall ∆ = 0 ⇔ the gros facteur and all its partial derivatives are 0

Voronoi diagrams of 3 lines

Proof: Main Lemma ∆ = 0 only if (y + a α) (a x + β) = 0

Recall ∆ = 0 ⇔ the gros facteur and all its partial derivatives are 0 Prove {∆ = 0, (y + a α) (a x + β) = 0} has no real solution

Voronoi diagrams of 3 lines

Proof: the gros facteur, all its partial derivatives, and Main Lemma ∆ = 0 only if (y + a α) (a x + β) = 0

Recall ∆ = 0 ⇔ the gros facteur and all its partial derivatives are 0

Consider the system of:

Prove {∆ = 0, (y + a α) (a x + β) = 0} has no real solution

Voronoi diagrams of 3 lines

Proof: the gros facteur, all its partial derivatives, and 1 − u (y + a α) = 0 1 − v (a x + β) = 0 Main Lemma ∆ = 0 only if (y + a α) (a x + β) = 0

Recall ∆ = 0 ⇔ the gros facteur and all its partial derivatives are 0

Consider the system of:

Prove {∆ = 0, (y + a α) (a x + β) = 0} has no real solution

Voronoi diagrams of 3 lines

Proof: the gros facteur, all its partial derivatives, and 1 − u (y + a α) = 0 1 − v (a x + β) = 0 1 − w (1 + α2 + β2) = 0 Main Lemma ∆ = 0 only if (y + a α) (a x + β) = 0

Recall ∆ = 0 ⇔ the gros facteur and all its partial derivatives are 0

Consider the system of:

Prove {∆ = 0, (y + a α) (a x + β) = 0} has no real solution

Voronoi diagrams of 3 lines

Proof: the gros facteur, all its partial derivatives, and 1 − u (y + a α) = 0 1 − v (a x + β) = 0 1 − w (1 + α2 + β2) = 0 1 − t Γ = 0

Γ = ` 2 a (yα − β x) − a2 + 1 ´2 + 3 (ax + β)2 + 3 a2 (y + aα)2 + 3 ` 1 + a2´2

Main Lemma ∆ = 0 only if (y + a α) (a x + β) = 0

Recall ∆ = 0 ⇔ the gros facteur and all its partial derivatives are 0

Consider the system of:

Prove {∆ = 0, (y + a α) (a x + β) = 0} has no real solution

Voronoi diagrams of 3 lines

Proof: the gros facteur, all its partial derivatives, and 1 − u (y + a α) = 0 1 − v (a x + β) = 0 1 − w (1 + α2 + β2) = 0 1 − t Γ = 0

Γ = ` 2 a (yα − β x) − a2 + 1 ´2 + 3 (ax + β)2 + 3 a2 (y + aα)2 + 3 ` 1 + a2´2

The Gr¨

• bner basis is 1, hence the system has no solution

(15mn with FGb [Faug` ere] in Maple)

Main Lemma ∆ = 0 only if (y + a α) (a x + β) = 0

Recall ∆ = 0 ⇔ the gros facteur and all its partial derivatives are 0

Consider the system of:

Prove {∆ = 0, (y + a α) (a x + β) = 0} has no real solution

Voronoi diagrams of 3 lines

About the proof. Start with the system of the gros facteur and all its partial derivatives

Voronoi diagrams of 3 lines

About the proof. If PQ is in basis B then study B ∪ {P} and B ∪ {1 − uP} Start with the system of the gros facteur and all its partial derivatives Compute reccursively Gr¨

• bner bases

Voronoi diagrams of 3 lines

About the proof. If PQ is in basis B then study B ∪ {P} and B ∪ {1 − uP} Set a = 2 Start with the system of the gros facteur and all its partial derivatives Compute reccursively Gr¨

• bner bases

Voronoi diagrams of 3 lines

About the proof. If PQ is in basis B then study B ∪ {P} and B ∪ {1 − uP} Set a = 2 By such manipulation, we add to the system Start with the system of the gros facteur and all its partial derivatives Compute reccursively Gr¨

• bner bases

1 − u (y + a α) = 0 1 − v (a x + β) = 0 1 − w (1 + α2 + β2) = 0

Voronoi diagrams of 3 lines

About the proof. If PQ is in basis B then study B ∪ {P} and B ∪ {1 − uP} Set a = 2 By such manipulation, we add to the system Start with the system of the gros facteur and all its partial derivatives Compute reccursively Gr¨

• bner bases

and eventually obtain in the basis

1 − u (y + a α) = 0 1 − v (a x + β) = 0 1 − w (1 + α2 + β2) = 0 Γ2 = (4 y α − 4 β x − 3)2 + 3 (2 x + β)2 + 12 (y + 2 α)2 + 75

Voronoi diagrams of 3 lines

About the proof. If PQ is in basis B then study B ∪ {P} and B ∪ {1 − uP} Set a = 2 By such manipulation, we add to the system Start with the system of the gros facteur and all its partial derivatives Compute reccursively Gr¨

• bner bases

and eventually obtain in the basis

1 − u (y + a α) = 0 1 − v (a x + β) = 0 1 − w (1 + α2 + β2) = 0 Γ2 = (4 y α − 4 β x − 3)2 + 3 (2 x + β)2 + 12 (y + 2 α)2 + 75

Repeat for a = 3, 5 and guess that

Γ = ` 2 a (yα − β x) − a2 + 1 ´2 + 3 (ax + β)2 + 3 a2 (y + aα)2 + 3 ` 1 + a2´2

Main Lemma ∆ = 0 only if (y + a α) = 0 or (a x + β) = 0

Voronoi diagrams of 3 lines

Main Lemma ∆ = 0 only if (y + a α) = 0 or (a x + β) = 0

Voronoi diagrams of 3 lines

When ∆ = 0, we prove that D(λ) > 0 thus D(λ) = 0 has 2 imaginary double roots λ1, λ2 We prove that Q1,2 + λiQ1,3 has rank 3 Classification of quadric intersections [DLLP08] ⇒ Trisector: cubic and line that do not intersect in P3(R)

(RAGLib) (FBb)

Main Lemma ∆ = 0 only if (y + a α) = 0 or (a x + β) = 0

Voronoi diagrams of 3 lines

When ∆ = 0, we prove that D(λ) > 0 thus D(λ) = 0 has 2 imaginary double roots λ1, λ2 We prove that Q1,2 + λiQ1,3 has rank 3 Classification of quadric intersections [DLLP08] ⇒ Trisector: cubic and line that do not intersect in P3(R) When ∆ > 0, the trisector is a smooth quartic The trisector has no singular point in P3(R)

(RAGLib) (FBb)

Voronoi diagrams of 3 lines

• The trisector has no singular point in P3(R)
• The trisector has no multiple points at infinity

❀ It is homeomorphic to 4 lines with no intersection point

• The set of triplets of lines in general position is connected
• The topology of the bisectors is invariant

The topology of the Voronoi diagram is invariant

Voronoi diagrams of 3 lines

The insight on the geometry yields a more direct proof gros facteur=0 ⇐ ⇒ the trisector contains a line

Voronoi diagrams of 3 lines

The insight on the geometry yields a more direct proof gros facteur=0 ⇐ ⇒ the trisector contains a line ⇐ ⇒ ℓ1, ℓ2, ℓ3 lie on a hyperboloid of revolution

Voronoi diagrams of 3 lines

The insight on the geometry yields a more direct proof gros facteur=0 ⇐ ⇒ the trisector contains a line ⇐ ⇒ ℓ1, ℓ2, ℓ3 lie on a hyperboloid of revolution ⇐ ⇒ the discriminant of the characteristic polynomial of the 3×3 top-left submatrix of the matrix representation

• f the hyperboloid containing ℓ1, ℓ2, ℓ3 is 0

The gros facteur is this discriminant

Plan

I.a. Properties of tangent lines in 3D

• II. Avoiding algebraic numbers

Intersections of two quadrics

• I. Computer algebra & theorems in geometry

I.b. Voronoi diagrams of three lines

Intersection of two quadrics

Challenge: Compute simple and usable exact parameterizations

• f the intersection of two quadrics

Perhaps one of the most basic problems in CAGD

From 1975 onwards: long list of publications but no satisfactory solution

Naive method

• Find a parameterization of Q1 in terms of u and v

Naive method

• Find a parameterization of Q1 in terms of u and v
• Insert the parameterization in the other quadric

❀ equation E in u and v

Naive method

• Find a parameterization of Q1 in terms of u and v
• Insert the parameterization in the other quadric

❀ equation E in u and v

• Solve the equation in v

Naive method

• Find a parameterization of Q1 in terms of u and v
• Insert the parameterization in the other quadric

❀ equation E in u and v

• Solve the equation in v
• Substitute v(u) in the parameterization

Naive method

• Find a parameterization of Q1 in terms of u and v
• Insert the parameterization in the other quadric

❀ equation E in u and v

• Solve the equation in v
• Substitute v(u) in the parameterization

Problem : Equation E may be of degree 4 in v (and u) not practical to handle

Past work

• Levin, A parametric algorithm for drawing pictures of solid
• bjects composed of quadric surfaces, 1976
• Algebraic approach : Sarraga (83),

Farouki, Neff & O’Connor (89), Wilf & Manor(93)

• Specific input/output : Shene & Johnstone (91, 92, 94)

Miller & Goldman (95)

• Other : Wang, Joe & Goldman (02)

Levin’s method (1976)

Idea: Intersection of two circles in the plane C1 : (x − a1)2 + (y − b1)2 − R2

1 = 0

C2 : (x − a2)2 + (y − b2)2 − R2

2 = 0

Levin’s method (1976)

Idea: Intersection of two circles in the plane C1 : (x − a1)2 + (y − b1)2 − R2

1 = 0

C2 : (x − a2)2 + (y − b2)2 − R2

2 = 0

Pencil: set of all linear combinations of equations C1 and C2 The intersection is independant of the choice of two distinct elements in the pencil

Levin’s method (1976)

Idea: Intersection of two circles in the plane C1 : (x − a1)2 + (y − b1)2 − R2

1 = 0

C2 : (x − a2)2 + (y − b2)2 − R2

2 = 0

Pencil: set of all linear combinations of equations C1 and C2 The intersection is independant of the choice of two distinct elements in the pencil Find a simple element in the pencil: the line C1 − C2

Levin’s method (1976)

Levin’s Main Theorem: In a pencil of quadrics, the simple quadrics are the ruled ones

(a quadric is ruled if it is a union of lines)

• There exists a ruled quadric in a pencil of quadrics
• Let λi be the solutions of the degree-three equation

det([P + λQ]3×3) = 0 One of the P + λiQ is ruled

Levin’s method (1976)

• find a ruled quadric in the pencil

Levin’s method (1976)

• find a ruled quadric in the pencil
• find a canonical form of this ruled quadric
• find a parameterization in u and v, linear in v

Levin’s method (1976)

• find a ruled quadric in the pencil
• find a canonical form of this ruled quadric
• find a parameterization in u and v, linear in v
• insert the parameterization in the other quadric

❀ equation in u and v,

• solve the equation of degree 2 in v

Levin’s method (1976)

• find a ruled quadric in the pencil
• find a canonical form of this ruled quadric
• find a parameterization in u and v, linear in v
• insert the parameterization in the other quadric

❀ equation in u and v,

• solve the equation of degree 2 in v
• substitute v(u) in the parameterization

Levin’s method (1976)

• find a ruled quadric in the pencil
• find a canonical form of this ruled quadric
• find a parameterization in u and v, linear in v
• insert the parameterization in the other quadric

❀ equation in u and v,

• solve the equation of degree 2 in v
• substitute v(u) in the parameterization

coeffs:

3

· + √ +

3

· + √

Levin’s method (1976)

• find a ruled quadric in the pencil
• find a canonical form of this ruled quadric
• find a parameterization in u and v, linear in v
• insert the parameterization in the other quadric

❀ equation in u and v,

• solve the equation of degree 2 in v
• substitute v(u) in the parameterization

coeffs:

3

· + √ +

3

· + √

coeffs: r · + q · +

3

p· + √ +

3

p· + √ +

3

p· + √ +

3

p· + √

Levin’s method (1976)

• find a ruled quadric in the pencil
• find a canonical form of this ruled quadric
• find a parameterization in u and v, linear in v
• insert the parameterization in the other quadric

❀ equation in u and v,

• solve the equation of degree 2 in v
• substitute v(u) in the parameterization

coeffs:

3

· + √ +

3

· + √

coeffs: r · + q · +

3

p· + √ +

3

p· + √ +

3

p· + √ +

3

p· + √

v u u t· + s · + r · + 3 q · + √ + 3 q · + √ + 3 q · + √ + 3 q · + √ + s · + r · + 3 q · + √ + · · ·

Levin’s method (1976)

Generically 5 nested radicals in the coefficients

Levin’s method (1976)

Generically 5 nested radicals in the coefficients What about in practice?

Levin’s method (1976)

Generically 5 nested radicals in the coefficients What about in practice? Hyperboloid of two sheets Ellipsoid x2 − y2 + z2 − xy − y + 1 = 0 2 x2 + y2 + 20 z2 − 4 yz + 4 y − 20 z − 20 = 0

Levin’s method (1976)

Generically 5 nested radicals in the coefficients What about in practice? Hyperboloid of two sheets Ellipsoid x2 − y2 + z2 − xy − y + 1 = 0 2 x2 + y2 + 20 z2 − 4 yz + 4 y − 20 z − 20 = 0 ∆ = −0.16 u4 − 1.0 u3 − 16.0 u2 − 0.6 u − 1.2 Approximate coefficients: x(u) = −0.44 u3−0.43 u2−1.4 u−0.15±(−2.9 10−6+2.3 10−7 u)

√ ∆ 1.2 u2+0.17 u+0.05

Exact coefficients: ❀ 20 megabytes Up to 200 megabytes even on simple instances!!!

Main ideas

Work in real projective space Avoid (at all cost) the appearance of algebraic numbers Classify the type of the intersection and design algorithms for every type of intersection

Our approach

Algorithm

• 1. Find in the pencil of QS and QT

a ruled quadric with rational coefficients

Algorithm

• 1. Find in the pencil of QS and QT

QS+λT is ruled ⇔ det(S + λT) 0 a ruled quadric with rational coefficients λ

det(S + λT)

λ

det(S + λT)

λ

det(S + λT)

λ0 ∈ Q λ0 ∈ Q λ0 ∈ Q

Algorithm

• 1. Find in the pencil of QS and QT

QS+λT is ruled ⇔ det(S + λT) 0 a ruled quadric with rational coefficients λ

det(S + λT)

λ

det(S + λT)

λ

det(S + λT)

λ0 ∈ Q λ0 ∈ Q λ0 ∈ Q

• Th. If ∀λ ∈ Q

det(S + λT) < 0 then QS ∩ QT = 2 points

• Proof. Based on canonical block diagonalization of

pairs of real symetric matrices [Uhlig (73)]

Algorithm

• 1. Find a ruled quadric QR with rational coefficients
• 2. Find a canonical form of this quadric

Algorithm

• 1. Find a ruled quadric QR with rational coefficients
• 2. Find a canonical form of this quadric

Use Gauss reduction of quadratic forms e.g. 2x2 + 4xy + 5y2 = 2(x + y)2 + 3y2 X = x + y Y = y for diagonalizing R

Algorithm

• 1. Find a ruled quadric QR with rational coefficients
• 2. Find a canonical form of this quadric

Use Gauss reduction of quadratic forms e.g. 2x2 + 4xy + 5y2 = 2(x + y)2 + 3y2 X = x + y Y = y for diagonalizing R To this point, everything is rational

Algorithm

• 1. Find a ruled quadric QR with rational coefficients
• 2. Find a canonical form of this quadric

Use Gauss reduction of quadratic forms e.g. 2x2 + 4xy + 5y2 = 2(x + y)2 + 3y2 X = x + y Y = y for diagonalizing R Focus on the generic case. det R > 0, thus we get

• r aX2 + bY 2 − cZ2 − dW 2 = 0

aX2 + bY 2 + cZ2 + dW 2 = 0 ⇒ QS ∩ QT = ∅ To this point, everything is rational

Algorithm

• 1. Find a ruled quadric QR with rational coefficients
• 2. Find a canonical form of this quadric
• 3. Parameterize aX2 + bY 2 − cZ2 − dW 2 = 0

❀ coeffs in Q(√ac, √ abcd)

(generic case)

[ ut+avs

a

, us−bvt

b

, ut−avs

√ac , us+bvt √ bd ], (u, v), (s, t) ∈ P1

Algorithm

• 1. Find a ruled quadric QR with rational coefficients
• 2. Find a canonical form of this quadric
• 3. Parameterize aX2 + bY 2 − cZ2 − dW 2 = 0

❀ coeffs in Q(√ac, √ abcd)

(generic case)

[ ut+avs

a

, us−bvt

b

, ut−avs

√ac , us+bvt √ bd ], (u, v), (s, t) ∈ P1

• Th. Q(

√ abcd) cannot be avoided

Algorithm

• 1. Find a ruled quadric QR with rational coefficients
• 2. Find a canonical form of this quadric
• 3. Parameterize aX2 + bY 2 − cZ2 − dW 2 = 0

❀ coeffs in Q(√ac, √ abcd)

(generic case)

[ ut+avs

a

, us−bvt

b

, ut−avs

√ac , us+bvt √ bd ], (u, v), (s, t) ∈ P1

• Th. Q(

√ abcd) cannot be avoided

• Th. Finding a parameterization over Q(

√ abcd)

• Finding a rational point on the quadric

Algorithm

• 1. Find a ruled quadric QR with rational coefficients
• 2. Find a canonical form of this quadric
• 3. Parameterize aX2 + bY 2 − cZ2 − dW 2 = 0

❀ coeffs in Q(√ac, √ abcd)

(generic case)

[ ut+avs

a

, us−bvt

b

, ut−avs

√ac , us+bvt √ bd ], (u, v), (s, t) ∈ P1

• Th. Q(

√ abcd) cannot be avoided

• Th. Finding a parameterization over Q(

√ abcd)

• Finding a rational point on the quadric

Param is worst-case optimal: ex. x2 + y2 − 3z2 − 11w2 = 0

Algorithm

• 1. Find a ruled quadric QR with rational coefficients

(generic case)

1’. Find a real point on QR Approximate it by a rational point p Consider the quadric QR′ through p

λ

det(S + λT)

Algorithm

• 1. Find a ruled quadric QR with rational coefficients
• 2. Find a canonical form of this quadric

(generic case)

1’. Find a real point on QR Approximate it by a rational point p Consider the quadric QR′ through p ❀ 4xy + z2 − ξw2 = 0

λ

det(S + λT)

Algorithm

• 1. Find a ruled quadric QR with rational coefficients
• 2. Find a canonical form of this quadric
• 3. Parameterize it ❀ coeffs in Q(√ξ)

(generic case)

1’. Find a real point on QR Approximate it by a rational point p Consider the quadric QR′ through p ❀ 4xy + z2 − ξw2 = 0

λ

det(S + λT)

Algorithm

• 1. Find a ruled quadric QR with rational coefficients
• 2. Find a canonical form of this quadric
• 3. Parameterize it ❀ coeffs in Q(√ξ)

(generic case)

1’. Find a real point on QR Approximate it by a rational point p Consider the quadric QR′ through p ❀ 4xy + z2 − ξw2 = 0

• 4. Substitute in the other quadric and solve

λ

det(S + λT)

Output

In the generic case (smooth quartic)

• X(u) = X1(u) ± X2(u)
• ∆(u) (optimal)
• X1, X2 and ∆ of degree 3, 1 and 4
• Coefficients in Q(√ξ), ξ ∈ N (worst-case optimal)

The approach in the generic case always works However

• parameterization may not be rational
• polynomial factorization with coeffs in Q(√ξ)

Specialized algorithms

• Non-generic intersection ⇐

⇒ det(S + λT) = 0 has multiple root(s)

• Characterize the real the type of the intersection by

the number and mutiplicity of multiple roots + inertia of associated matrices ❀ Use that geometric information ❀ Use the quadric QµS+λT associated to the mult. root to write dedicated algorithms to parameterize the intersection

Theory Implementation

QI: efficient C++ implementation

(∼ 20 000 lines)

On-line web server

Intersection of two quadrics

[Dupont, Lazard, Lazard, Petitjean, J. of Symb. Comput. 2008] [Lazard, Pe˜ naranda, Petitjean, Comput. Geom. Theory & App., 2006]

Classification of the type of intersection in P3(R)

(smooth quartic, cubic & line, 2 conics, etc.)

Algorithm for parameterizing each component Rational parameterization if one exists Small coefficients (near-optimal degree of the extension of Q)

Experiments

Random quadrics, coeffs up to 10 digits: < 50 ms Chess set: 6 diff. pieces (108 quadrics), 971 intersetions 3.4 ms on average

Intersection of two quadrics

Implementation

C++ built on GMP and LiDIA

Intersection of two quadrics

An example: 4 x2 + z2 − w2 = 0 x2 + 4 y2 − z2 − w2 = 0

[Wang, Joe, Goldman, 2002]

Intersection of two quadrics

An example: 4 x2 + z2 − w2 = 0 x2 + 4 y2 − z2 − w2 = 0

X(u)= B B B B B B @ 0.0 1131.3708 u3−5760.0 u2+10861.1602 u−8192.0 −1600.0 u3+10861.1602 u2−21504.0 u+11585.2375 1600.0 u3+3620.2867 u2+5120.0 u+11585.2375 1 C C C C C C A ± B B B B B B @ −80.0 u+1181.0193 0.0 0.0 0.0 1 C C C C C C A

√

905.0967 u3−3328.0 u2+2896.3094 u

[Wang, Joe, Goldman, 2002]

(approximated)

Intersection of two quadrics

An example: 4 x2 + z2 − w2 = 0 x2 + 4 y2 − z2 − w2 = 0

X(u) = B B @ 2 u3 − 6 u 7 u2 + 3 10 u2 − 6 2 u3 + 18 u 1 C C A ± B B @ −2 u 2 u 2 1 C C A p −3 u4 + 26 u2 − 3. Running time < 10ms

Our approach (exact)

Conclusion

Algebraic computations in geometric algorithms Avoiding the appearance of algebraic numbers Using computer algebra tools in geometric algorithms Major progress can be made on specific instances Possible and efficient. . . still it’s not quite yet plug-and-play Advances in computer algebra ❀ possible to solve non-trivial geometric problems Computer algebra & theorems in geometry ❀ Theorems in geometry Still nothing is easy. . .