Pedro M M De Castro

Frédéric Cazals

Sébastien Loriot

Monique Teillaud

Design of the cgal 3D Spherical Kernel and application to arrangements of circles on a sphere *

This paper presents a cgal kernel for algorithms manipulating 3D spheres, circles, and circular arcs. The paper makes three contributions. First, the mathematics underlying two non trivial predicates are presented. Second, the design of the kernel concept is developed, and the connexion between the mathematics and this design is established. In particular, we show how two different frameworks can be combined: one for the general setting, and one dedicated to the case where all the objects handled lie on a reference sphere. Finally, an assessment about the efficacy of the 3D Spherical Kernel is made through the calculation of the exact arrangement of circles on a sphere. On average while computing arrangements with few degeneracies (on sample molecular models), it is shown that certifying the result incurs a modest factor of two with respect to calculations using a plain double arithmetic.

Introduction

1.1 Linear versus Curved Primitives in Computational Geometry.

Until recently geometric algorithms used to limit their framework to linear objects in affine Euclidean space (points, segments, triangles, . . . ), curved objects being discretized by linear elements. However, handling directly curved objects allows one to skip this discretization process, thus reducing the number of primitives and yielding algorithms with better combinatorial complexity. For these reasons, the direct manipulation of curved objects is an important challenge in computational geometry. Of particular interest in the realm of curved objects are the simplest primitives, namely circles, circular arcs and spheres, which are ubiquitous to either represent or approximate more elaborate objects.

Regarding one dimensional primitives, circular arcs in the plane are involved in the design and the manufacturing of printed and integrated circuits [START_REF] De Castro | Exact and efficient computations on circles in CGAL[END_REF], and they also generate an increasing interest for approximating curves [START_REF] Drysdale | Approximation of an open polygonal curve with a minimum number of circular arcs and biarcs[END_REF][START_REF] Aichholzer | Computational and structural advantages of circular boundary representation[END_REF]. Balls and spheres are probably even more central, especially to approximate 3D objects. Given a known 3D object, balls can be used to provide hierarchical approximations, for example to efficiently perform collision detection [START_REF] Agarwal | Collision detection for deforming necklaces[END_REF] or multi-scale visualization [START_REF] Rusinkiewicz | QSplat: a multiresolution point rendering system for large meshes[END_REF]. Given a 3D point cloud sampling an unknown object, balls are also key for the inference of geometric and topological properties of the unknown object from the samples. Topologically speaking, balls can indeed be used to reconstruct the object by mimicking the medial axis transform, which is the bottom-line of say the power crust algorithm [START_REF] Amenta | The power crust, unions of balls, and the medial axis transform[END_REF]. From a more geometric standpoint, a number of quantities can be inferred from balls centered at the sample points. One may mention the so-called boundary measure of a point cloud, from which singular points and sharp features of the sampled model can be estimated [START_REF] Chazal | Stability of boundary measure[END_REF]. This inference requires computing the volume of a ball contained within its Voronoi region, a calculation involving geometric constructions on a sphere [START_REF] Mcconkey | Quantification of protein surfaces, volumes and atomatom contacts using a constrained Voronoi procedure[END_REF][START_REF] Avis | Computing the volume of the union of spheres[END_REF]. Given a 3D polyhedron, additional geometric modeling constructions involving spheres are the computation of offset surfaces. (Notice that the offset surface calculation covers the 3D Minkowski sum between a polyhedron and a sphere.) Interested readers may also consult chapter 11.12 of the Visionbib bibliography project at http://www.visionbib.com/bibliography/describe482.html, which contains an exhaustive list of representations involving spheres.

Finally, spheres are also central in molecular modeling since a Van der Waals model features a collection of balls whose radii depend upon the atom type and its covalent environment. In particular, the combinatorial and geometric properties of the three main definitions of surface of a molecule used (the Van der Waals surface, the solvent-accessible surface and the molecular surface) read from the boundary of a union of balls [START_REF] Connolly | Molecular surfaces: a review[END_REF].

From Algorithms to Implementations

Predicates versus Constructions in the Exact Geometric Computation Paradigm. Recent implementations of geometric algorithms usually distinguish predicates and constructions. In standard folklore, a predicate is a test function returning an element of a discrete set of values (for instance: do the two curves c and c ′ intersect?). A construction is a function that constructs new objects (for instance: compute the intersection points between curves c and c ′ ). Since rounding errors due to floating point calculations often produces inconsistencies [START_REF] Kettner | Classroom examples of robustness problems in geometric computations[END_REF], the robustness of geometric algorithms relies on the exact evaluation of predicates [START_REF] Yap | The exact computation paradigm[END_REF]. This problem usually boils down to evaluating the sign of polynomial expressions on the input data, which can be done efficiently resorting to arithmetic filtering techniques [START_REF] Li | Recent progress in exact geometric computation[END_REF]. Aside this classical strategy, one further refinement consists of storing, in the course of the predicate evaluation, intermediate geometric objects to be reused later on. In the realm of curved objects, the former strategy typically consists of using an implicit encoding of an algebraic number by some polynomial vanishing on it, which calls for resultant-like calculations on the input polynomials. The latter relies on exact representation of algebraic numbers [START_REF]Core Number Library[END_REF][START_REF]Library for Efficient Data Types and Algorithms[END_REF], and requires in the end running multi-precision calculations down to separation bounds in the worst cases [START_REF] Burnikel | A separation bound for real algebraic expressions[END_REF][START_REF] Li | A new constructive root bound for algebraic expressions[END_REF]. But in doing so, one reuses intermediate calculations, and the design of predicates, a difficult task for algebraic objects, is made easier. In fact, work on filtered constructions [START_REF] Funke | LOOK: A Lazy Object-Oriented Kernel for geometric computation[END_REF][START_REF] Fabri | A generic lazy evaluation scheme for exact geometric computations[END_REF] may reconcile the evaluation of predicates and constructions. The implementation of such a framework still needs to be worked out to possibly surpass the approach relying on predicates without intermediate constructions.

The CGAL Library, Kernels and Traits Classes. The goal of the cgal open source project is to promote research in computational geometry, so as to make reference algorithms available as robust programs for academic research and industrial applications. We refer the reader to the cgal web site and to the literature for more background on the history and running of the cgal project [START_REF] Cgal | [END_REF][START_REF] Fogel | Generic programming and the CGAL library[END_REF]. Geometric algorithms are generic, and the geometric classes are instantiated with a Traits class-offering the minimum set of functionalities required. Such a template parameter is supposed to provide a set of functionalities, with prescribed signatures, described and documented as a concept. In the cgal library, constant size geometric primitive objects, together with general purpose predicates and constructions on them, are made available in kernels. As kernels can be directly used as traits classes by several cgal classes, cgal kernels are documented as C++ concepts-we say that one implementation of a given concept is a model of the concept.

Circles and Spheres: from Primitives to Algorithms

As a follow-up to algorithms geared towards linear primitives, a recent trend in computational geometry has been the design of algorithms and programs geared towards curves and surfaces. The work mostly focused on computing arrangements in 2D and 3D, which involves primitive computations such as intersection and relative position of objects.

As far as we know, in the planar case mapc was the first geometric library handling general algebraic curves [START_REF] Keyser | MAPC: a library for efficient and exact manipulation of algebraic points and curves[END_REF], but all degenerate cases were not handled. The Exacus prototype library allows to exactly compute arrangements of conics and cubics in the plane [START_REF] Exacus | Efficient and Exact Algorithms for Curves and Surfaces[END_REF]. The first version of a kernel for non-linear objects (circles and circular arcs in 2D) was released in cgal 3.2 [START_REF] Pion | 2D circular kernel[END_REF]. The cgal arrangement package can use this kernel to compute an arrangement of circular arcs in the plane; it can also deal with conics or Bézier curves [START_REF] Wein | 2D arrangements[END_REF].

In the three dimensional setting, Andrade and Stolfi proposed selected exact manipulations of 3D circular objects lying on a given sphere, with an implementation in Modula-3 [START_REF] Andrade | Exact algorithms for circles on the sphere[END_REF]. Their representation of objects uses Plücker coordinates and is restricted to the case of circular arcs on a given sphere, whereas our framework handles general circles and circular arcs in 3D. To the best of our knowledge, there exist two algorithms computing the exact arrangement induced by a set of circles on a sphere. The first one is generic as its specification handles general curves on a parametric surface, the strategy consisting of sweeping the parametric domain of the surface [START_REF] Berberich | Sweeping and maintaining twodimensional arrangements on surfaces: A first step[END_REF]. But for a surface which is a sphere, only great circles are supported in the current implementation. The second one is an algorithm which consists of sweeping a sphere with a meridian anchored at the poles. The algorithm handles all types of circles and any degenerate case [START_REF] Cazals | Computing the exact arrangement of circles on a sphere, with applications in structural biology[END_REF]. We report in section 5 how the 3D Spherical Kernel is used in the implementation of this algorithm.

While discussing arrangements and for the sake of completeness, one should also mention papers dealing with quadrics in 3D. The first complete and exact implementation computing a planar map induced by the intersection curves of a set of quadrics running on the surface of one of them, say Q, is described in [START_REF] Berberich | An exact, complete and efficient implementation for computing planar maps of quadric intersection curves[END_REF]. The algorithm reduces the problem to a planar arrangement of algebraic curves using a projection into a plane. Moreover, it is shown how this planar arrangement can be used to compute the arrangement on Q-a lifting step which can be avoided when directly computing the arrangement on Q in 3D. Finally, an algorithm to compute the exact volumetric decomposition induced by an arrangement of quadrics is developed in [START_REF] Mourrain | On the computation of an arrangement of quadrics in 3d[END_REF]. The strategy does not resort to surface arrangements, but no implementation is provided.

Contributions and Paper Overview

This work provides the first complete and efficient implementation of essential basic functionalities to manipulate in an exact way spheres, circles and circular arcs in 3D, so as to match most of the needs mentioned in section 1.1. As opposed to the general algorithms discussed in 1.3, the perspective is that of a hand-crafted package geared towards the aforementioned objects. Following best practices, a clear distinction is made between the algebraic and the geometric aspects on the one hand, and on the concepts of a kernel and its implementation on the other hand. The 3D Spherical Kernel concept is accompanied by an implementation, which is used to compute the exact arrangements of circles on a sphere. As a quality label, the basic part of the package, related with general spheres, circles and circular arcs in 3D, has been reviewed and accepted by the cgal Editorial Board [START_REF] De Castro | 3D spherical kernel[END_REF]. The second part of the package, offering functionalities on a reference sphere, will follow the process in a second step. In particular, the diffusion of our code will enable a number of endeavors. As an example, it will foster the development of a traits class providing the predicates required to handle all types of circles in the general sweep-line algorithm for curves on a surface [START_REF] Berberich | Sweeping and maintaining twodimensional arrangements on surfaces: A first step[END_REF].

The paper is organized as follows. Section 2 presents the mathematics underlying two non trivial predicates. Section 3 presents the design of the 3D Spherical Kernel and of the Algebraic Kernel. In Section 4, we bridge the gap between the geometric calculations presented in Section 2 and the functionalities of the kernel. The application to computing an arrangement of circles on a sphere is presented in Section 5. Finally, experiments on random arrangements and molecular models are reported in Section 6. Due to space limitations, the reader is referred to [START_REF] De Castro | Design of the CGAL spherical kernel and application to arrangements of circles on a sphere[END_REF] for additional details skipped in this paper -calculations for predicates and constructions as well as example C++ code.

Mathematical Background

This section focuses on general basics underlying the 3D Spherical Kernel, as well as on the mathematics needed for two non-trivial operations: The construction of the θ-extremal points of a circle on a given sphere, and the computation of the relative orientation of the tangents of two intersecting circular arcs on a given sphere.

Preliminaries and Notations

The center of a sphere S i is denoted c i = (x i , y i , z i ), and its radius r i . We assume the Cartesian coordinates of the center and the squared radius to be rational numbers. All objects (planes, lines, spheres, circles) whose equations have rational coefficients are termed rational.

The x-coordinate of point p and vector u are denoted p x and u x , and similarly for y-and z-coordinates. The dot and vector products of two vectors u and v are respectively denoted < u, v > and u ∧ v. The squared norm of vector u is denoted u 2 =< u, u >, and its norm u . The sign of a real number, denoted Sign(x), is such that Sign(x) ∈ {-1, 0, 1}. The vector defined by two points p and q is denoted pq. If o stands for the origin, and p is a point, the vector op is denoted p.

The power of point p w.r.t. sphere S i is defined by π(p, S i ) = pc 2 i -r 2 i . The radical plane RP ij of any two spheres S i and S j is the plane consisting of the points having equal power with respect to the two spheres. Its equation is 2 < p, c i c j > +c i 2 -c j 2 + r 2 j -r 2 i = 0. Whenever these two spheres intersect, their intersection circle C ij is also defined as the intersection between either sphere and their radical plane.

The center and radius of C ij are denoted c ij and r ij . It can be checked that

c i c ij = cic 2 j -r 2 j +r 2 i 2cic 2 j c i c j and r 2 ij = r 2 i -c i c 2 ij .
A root of a degree n polynomial with rational coefficients is called an algebraic number of degree n. The following observation shows that algebraic numbers of degree two play a central role in our predicates and constructions: Observation 1 The Cartesian coordinates of the intersection points of a rational sphere and a rational line are algebraic numbers of degree two in the same extension.

We also recall the following properties that will be used in our computations:

Observation 2 Consider two algebraic numbers a and b, of degree at most two in the same algebraic extension. Then a + b, a × b and 1/a (if a = 0) are algebraic numbers of degree at most two, and they also belong to the same algebraic extension.

Objects on a Reference Sphere

In the following and without loss of generality, we will consider a sphere S 0 centered at the origin called the reference sphere. S 0 is equipped with cylindrical coordinates (θ, z). Let M θ be a parametrized meridian on S 0 , with θ ∈ (0, 2π]. We first define a classification of circles on S 0 , and related quantities:

Definition 1 Consider a circle on a reference sphere S 0 . Such a circle is said to be polar (respectively bipolar) if it goes through one pole (respectively the two poles) of S 0 . A circle which separates S 0 into two connected components, each containing one pole, is termed threaded. Any other circle is termed normal. For polar and normal circles, we introduce the notion of θ-extremal points.

Definition 2

The θ-extremal points of a normal circle are the two points where a meridian of S 0 is tangent to the circle. The θ-extremal point of a polar circle is the pole the circle goes through.

For all but threaded circles, we define the notion of θ-extremal value(s): Definition 3 For a normal circle, the θ-extremal values are the θ-coordinates of its θ-extremal points. For a polar circle, the θ-extremal values are the two values of θ such that M θ is tangent to the circle at the pole. For a bipolar circle, the θ-extremal values are the two values of θ such that M θ is contained in the plane of the circle.

Notice that in standard analysis and geometric folklore, the θ-extremal points are also called critical points. Notice also that the θ-extremal points should not be confused with the endpoints of an arbitrary circular arc.

Definition 4 A circular arc on a reference sphere is said to be θ-monotone if its intersection with any half-plane bounded by the line going through the two poles is empty or reduces to one point.

Consider the four types of circles, as illustrated on Fig. 1. A circular arc on a normal circle is θ-monotone if it contains none of the θ-extremal points of its supporting circle, excepted maybe as an endpoint. The same remark is valid on a polar circle, considering its θ-extremal point. Any circular arc constructed on a threaded circle is by definition θ-monotone while no θ-monotone circular arc can be found on a bipolar circle.

If C 0i is a normal circle, the two circular arcs delimited by the θ-extremal points are called the upper or lower circular arcs as evidenced by their relative position along the z-axis. By extension, a θ-monotone circular arc included into the upper (lower) arc is also termed upper (lower) and a θ-monotone circular arc on a polar circle passing through the north (south) pole is termed lower (upper).

Computing the θ-extremal Points of a Normal Circle

We start with a proposition concerned with the z-coordinate of θ-extremal points:

Proposition 1 Consider a normal circle C 0i defined by the intersection of two rational spheres. The z-coordinate of its two θ-extremal points is the following rational number:

z = 2z i r 2 0 c i 2 + r 2 0 -r 2 i .
Proof. Recall that the reference sphere S 0 is centered at the origin, and let C 0i be a normal circle intersection of S 0 and S i . If the z-coordinate c 0iz of its center c 0i is equal to 0, a symmetry argument with respect to the plane z = 0 imposes that the z-coordinate of the θ-extremal points is also null. In the following, we therefore assume that c 0iz = 0. Consider the zero-dimensional intersection set between the circle C 0i (or equivalently the intersection S 0 ∩ RP 0i ) and the half-plane P (θ) defining the meridian M θ (i.e. M θ = S 0 ∩ P (θ)). The θ-extremal points of the circle correspond to the case where this intersection set reduces to one point. Denoting p = (x, y, z) such an intersection point, the corresponding polynomial system reads as:

     P (θ) : -x sin θ + y cos θ = 0 RP 0i : 2 < p, c 0 c i > +c 0 2 -c i 2 + r 2 i -r 2 0 = 0 S 0 : x 2 + y 2 + z 2 -r 2 0 = 0 (1)
Assuming x = 0, or equivalently θ = π 2 [π], we investigate this system using tan θ. (If x = 0, we work with cot θ by swapping the roles of x and y. Notice that x and y cannot vanish simultaneously for a normal circle, as such a circle does not contain a pole.) The previous system is tantamount to:

     y = x tan θ z = ci 2 +r 2 0 -r 2 i -2x(xi+yi tan θ) 2zi x 2 + x 2 tan 2 θ + (ci 2 +r 2 0 -r 2 i -2x(xi+yi tan θ)) 2 4z 2 i -r 2 0 = 0 (2)
Rewrite the last equation of the previous system

Ax 2 + Bx + C = 0, (3) with
       A = (1 + tan 2 θ) + (xi+yi tan θ) 2 z 2 i B = - (xi+yi tan θ)(ci 2 +r 2 0 -r 2 i ) z 2 i C = (ci 2 +r 2 0 -r 2 i ) 2 4z 2 i -r 2 0 (4)
The values of θ sought are such that Eq. ( 3) has a single solution, namely x = -B/(2A). Imposing that the discriminant of this polynomial vanishes yields the condition D = 0, with D = B 2 -4AC, whence tan θ as an algebraic number of degree two. Letting T stand for tan θ and denoting λ i = c i 2 + r 2 0 -r 2 i , we have:

D = D2 (-λ 2 i + 4r 2 0 z 2 i + 4r 2 0 y 2 i ) T 2 + D1 (8x i y i r 2 0 ) T + D0 (-λ 2 i + 4r 2 0 z 2 i + 4r 2 0 x 2 i ) (5)
Let us investigate the cases where polynomial D has degree two and one. As we shall see, in the first case, no θ-extremal point lies in the plane x = 0, while the second one corresponds to the situation where one θ-extremal point lies in plane x = 0. The latter case is left to the reader and we describe here only the former one. We assume that D 2 = 0: Polynomial D is of degree two. We first make the connexion between the discriminant δ of D and the geometry, and proceed with the value of the z-coordinate of θ-extremal points. The discriminant δ of D may be written as:

δ = 4 δ1 (-2r 0 z i + λ i ) δ2 (2r 0 z i + λ i ) δ3 (4r 2 0 c i 2 -λ i ) (6)
Since the power of a pole of S 0 w.r.t. to sphere S i reads as

(c i -(0, 0, ±r 0 )) 2 -r 2 i = c i 2 + r 2 0 -r 2 i ± 2z i r 0 = ±2r 0 z i + λ i , polynomial δ 1 (respectively δ 2 )
is the power of the north (respectively south) pole w.r.t. S i . The power of a point w.r.t. to a sphere is negative (respectively positive, zero) iff the point is inside (respectively outside, on) the sphere. Therefore, δ 1 is negative (respectively positive, zero) if the north pole is inside (respectively outside, on) S i . The same observation holds for δ 2 and the south pole. Finally, for δ 3 , observe that:

δ 3 = δ3 1 (-c i 2 + (r 0 + r i ) 2 ) δ3 2 (c i 2 -(r 0 -r i ) 2 ) . (7)
The combinations of signs of δ 31 and δ 32 are reported in Table 1, and the discussion of the sign of δ is summarized in Table 2. These tables prove that as expected, one has δ > 0 for all normal circles, so that Eq. ( 5) can be cancelled in two different ways corresponding to the two θ-extremal points. Finally, substituting any solution of Eq. ( 5) into the expression of z from the system (2) proves the claim under the assumption D 2 = 0. Notice that c i 2 + r 2 0 -r 2 i = 0 corresponds to the case where the circle is bipolar. 7). Abusing terminology, S i inside S j stands for S i inside the ball associated to S j while S i outside S j stands for the balls associated to S i and S j are disjoint.

δ 1 δ 2 δ 3 δ circle type + + + + normal circle + + - - S 0 and S i do not intersect + - + - threaded circle - + + - threaded circle + - - + impossible - + - + impossible - - + + normal circle - - - - S 0 and S i do not intersect 0 = 0 ≥ 0 0 polar circle if δ 3 > 0, S 0 and S i are tangent otherwise = 0 0 ≥ 0 0 polar circle if δ 3 > 0, S 0 and S i are tangent otherwise 0 = 0 - 0 impossible = 0 0 - 0 impossible ± ± 0 0 tangent spheres if δ 1 δ 2 > 0, impossible otherwise 0 0 {±, 0} 0 bipolar circle if δ 3 > 0, impossible otherwise
Table 2: Sign of δ from Eq. ( 6). δ 1 (respectively δ 2 ) is the power of the north (respectively south) pole w.r.t. S i while δ 3 indicates whether S 0 intersects S i .

Proposition 1 has an interesting consequence:

Corollary 1 The θ-extremal points of a normal circle C 0i are defined by the intersection between the reference sphere, and the intersection line between the radical plane RP 0i and the horizontal plane defined by z =

2zir 2 0 ci 2 +r 2 0 -r 2 i
. Therefore, by Observation 1, the x-and y-coordinates of these points are algebraic numbers of degree two in the same extension.

Computing the Relative Orientation of the Tangents of two Circles

Consider two circles, each supported by a rational plane, and intersecting transversally at point p on S 0 . Under the assumptions that none of the circles is bipolar and that p is neither a pole nor a θ-extremal point of a normal circle, we wish to compute the relative orientation of the tangents of these two circles at p.

Under the previous assumptions, denote t k , k = i or j, the tangent vector associated to circle C 0k at point p, as illustrated on Fig. 2. Tangents vectors t i and t j are chosen so as to point towards increasing values of θ, locally in the tangent space of S 0 at p. Because t i and t j are orthogonal to p, the relative orientation of the tangents is given by Sign(∆), with

∆ =< t i ∧ t j , p > . (8)
The details are as follows. The tangent to a circle is supported by the line intersection of the plane of the circle together with the tangent plane of S 0 at p. More precisely, consider the vectors n i = c 0i p and n j = c 0j p. For k = i, j the tangent vectors are obtained as follows: For a normal or a polar circle, t k = β k c 0k ∧ n k , with β k = -1 (respectively 1) if the circular arc containing p is an upper (respectively a lower) one; For a threaded circle, let m k be c 0k if the circle is not a great circle and otherwise any normal vector of the plane of the circle, then

t k = γ k (m k ∧ n k ), with γ k = ±1 such that γ k m kz > 0. S 0 c 0i c 0j C 0i C 0j n i n j -t i t j p z M θ o n
Figure 2: Tangent vectors t i and t j to circles C 0i and C 0j at one of their intersection points p. Vector n is the normal to S 0 at p.

The Cartesian coordinates of p are three algebraic numbers of degree two in the same extension. Computing tangent vectors using the formulation t k = ±(c 0ky p z -c 0kz p y , c 0kz p x -c 0kx p z , c 0kx p yc 0ky p x ) t , we now give the cost of computing ∆ in term of products Π T ;T ′ and sums Σ T ;T ′ , where T and T ′ are either algebraic numbers of degree two (denoted AN op ) or rational numbers (denoted R). AN op is a trivariate rational polynomial expression in the coordinates of p. This ensures that it is an algebraic number of degree two lying in the same extension as p. If considering the tree of arithmetic operations representing such an expression, the subscript op is equal to the maximum number of edges between the final expression and a coordinate of p. For instance each coordinate of t k is a AN 2 obtained through two Π R;AN0 followed by one Σ AN1;AN1 . A naive computation of ∆ requires two Π AN2;AN2 and one Σ AN3;AN3 for each coordinate of the vector product; three Π AN4;AN0 , one Σ AN5;AN5 and one Σ AN6;AN5 for the dot product. In Observation 4, omitting products of rationals, we show that if no circle is a great circle, we can actually compute the sign of ∆ with only two Π AN2;AN2 , two Π R;AN3 and one Σ AN4;AN4 . We use the following elementary observation:

Observation 3 Let ∆ ′ = tj z tj -ti z
ti . The sign of ∆ matches that of ∆ ′ .

Notice that the sign of ∆ ′ is trivially inferred in the following three cases: t iz = 0 and t j z = 0; t iz = 0 or t j z = 0; t iz and t j z have different signs.

Otherwise we have:

Observation 4 If none of the circles involved is a great circle and the signs of t iz and t j z are identical, the sign of ∆ ′ is given by:

Sign(∆ ′ ) = Sign(t j z )Sign c 0i 2 r 2 0i (t j z ) 2 -c 0j 2 r 2 0j (t iz ) 2
Proof. In the following, we suppose without loss of generality that the signs of t iz and t j z are identical and non null. The sign of ∆ ′ is the same as that of t i t j z -t j t iz . and we have:

t 2 i t j 2 z -t 2 j t i 2 z = ( t i t j z + t j t iz )( t i t j z -t j t iz ) (9)
And since we supposed that Sign(t j z ) = Sign(t iz ), we also have

Sign(∆ ′ ) = Sign(t j z )Sign(t 2 i (t j z ) 2 -t 2 j (t iz ) 2 ). ( 10
)
Recall that the Cartesian coordinates of the center of circle C 0k , k = i, j, are rational numbers. As c 0k and c 0k p are orthogonal, the sine of the angle between the vectors is 1. Since t k = ±c 0k ∧ c 0k p, we have t 2 k = c 0k 2 c 0k p 2 , which we can rewrite t 2 k = c 0k 2 r 2 0k .

Software Design

In this section, we sketch the design perspective of cgal kernels in general, and focus on the 3D Spherical Kernel.

CGAL Kernels: Rationale

The general design of a kernel concept for curved objects initially proposed in [START_REF] Emiris | Towards an open curved kernel[END_REF] is driven by the following goals: (i)interoperability of any model of the kernel concept with cgal geometric algorithms;

(ii)re-usability of the existing cgal kernel models for linear objects; (iii)genericity and flexibility: ability to use other linear kernels than the cgal ones, and independence from a particular implementation of the algebraic operations needed. The declaration of a kernel for curved objects is the following:

template < typename LinearKernel, typename AlgebraicKernel > class Curved_kernel;

Practically, to indicate whether an object is a member of the 3D Spherical Kernel or of the Algebraic Kernel, we use the following prefixes SK:: and AK:: respectively.

Let us now comment briefly on the two template arguments. Using the extensibility and adaptability scheme of the cgal kernel [START_REF] Hert | An adaptable and extensible geometry kernel[END_REF], the Curved kernel inherits from the model of LinearKernel given as template parameter, which allows it to directly benefit from all the functionalities on linear objects. The LinearKernel concept will not be presented here. It coincides with the basic cgal kernel concept extensively documented in the cgal manual. The requirements on the AlgebraicKernel concept listed in Section 3.3 are guided by the functionalities offered by the 3D Spherical Kernel (Section 3.2). Notice that this Algebraic Kernel concept is actually called AlgebraicKernelForSpheres in the cgal package because it is restricted to functionalities required by the 3D Spherical Kernel. Nevertheless, we use the shorter name AlgebraicKernel in this paper.

At the embedding level, the number type used to represent the objects of the linear kernel (coordinates for points, coefficients of equations for planes, spheres) is supposed to lie in a field number type (a type providing elementary operations +, -, ×, /), that will typically be the rationals. From now on and for the sake of simplicity, we shall refer to this basic number type as rational.

Regarding algebra, the key idea in this design is that each primitive object of the 3D Spherical Kernel is internally represented by a corresponding object of the Algebraic Kernel-its equation or coordinates. Moreover, each functionality of the 3D Spherical Kernel is obtained by a combination of calls of functionalities of the Algebraic Kernel on the algebraic representations of its arguments. Fig. 3 illustrates how algebra and geometry communicate, on the example of the computation of the intersection of a circle and a plane: We first retrieve the underlying primitive geometric objects, get their equations as polynomials of the Algebraic Kernel, and solve the system of corresponding equations using the solver provided by the Algebraic Kernel. Finally, the root of the system is used to retrieve the geometric intersection.

3D Spherical Kernel

To list the main geometric functionalities offered upon instantiation of a model of the 3D Spherical Kernel concept, we discuss in turn the types provided, the predicates and constructions, and the connexion between geometry and algebra.

Types. As seen from Table 3, the types decompose into three groups. The first four types are inherited from the 3D Linear Kernel, whereas the other objects are introduced by the 3D Spherical Kernel itself. Then, the 3D Spherical Kernel concept can be decomposed into two different parts: The first one considers general objects, while the second part introduces specific concepts for objects lying on a given reference sphere-as evidenced by the suffix ORS which stands for OnReferenceSphere. Note that the concept ThetaRep, designed to represent the θ-coordinate of a point, requires no accessor. It must only be comparable through function object CompareTheta 3. Note that a circular arc can be defined by a supporting circle and two endpoints. It is unambiguously defined as the set of points lying on the circle, when walking counterclockwise from its source to its target in the positive plane containing the circle-a plane is positive if its equation is of the form ax+by+cz+d = 0 with (a, b, c) > (0, 0, 0) according to the lexicographic order. Predicates and Constructions. Predicates and constructions are provided by the kernel as C++ functions objects (called functors). Note that a functor can be used as a template argument.

(A m i , A m j , M) Order along M {<,=,>} (p r , A m i ) p r vs. A m i along M θ ∋ p r SK::CompareZToRight 3 {<,=,>} (A m i , A m j )
Right of common pt order Notations: p, A, P , M, S, θ, w, stand respectively for a point, a circular arc, a plane, a meridian included in a rational plane, a sphere, a θ-value of type ThetaRep, any object of the 3D Spherical Kernel. Two instances of any type are distinguished using a subscript i and j. A superscript m on a circular arc indicates that this arc is θ-monotone-obviously on a reference sphere. A superscript r on a point indicates that such a point is on a reference sphere. Notice that for equality test, w i and w j are of the same type.

Table 4: Predicates of the 3D Spherical Kernel.

A model of the 3D Spherical Kernel must provide the functionality described in the sequel. We list only the most important requirements in Table 4 (predicates) and Table 5 (constructions). A construction returns an iterator over a (possibly empty) set of solutions associated to the problem solved by the function. Communication with Algebra. Functor SK::GetEquation returns a polynomial of the Algebraic Kernel providing an equation of a plane or a sphere. As explained at the beginning of Section 3, this is the essential link between the 3D Spherical Kernel and the Algebraic Kernel.

Algebraic Kernel for Spheres

The operations provided by the cgal 3D Spherical Kernel make heavy use of algebraic operations. The AlgebraicKernel parameter has a crucial role in particular for the robustness of the 3D Spherical Kernel.

Construction

Arguments Returned objects SK::Intersect 3

Two or three primitives Their intersection(s) SK::ThetaExtremePoint 3 A circle or a circular arc, Its θ-extremal point(s) on a reference sphere SK::MakeThetaMonotone 3 A circle or a circular arc, Its θ-monotone on a reference sphere sub-arc(s)

Table 5: Constructions of the 3D Spherical Kernel.

Sorting Circular Arcs at a Common Point. Let p be an intersection point of two θ-monotone circular arcs and suppose that p is on a meridian M θ . The predicate SK::CompareZToRight 3 consists of finding the relative position (above, on or below) of the two arcs to the right of p, i.e. for the angle value θ + ε with ε arbitrarily small. Notice that the predicate is not defined if point p is a pole, and that the associated circles are of any type but bipolar as no θ-monotone circular arc is defined on such circles. The two supporting circles are denoted C 0i and C 0j . We now consider three exclusive cases:

-If p matches a θ-extremal point of at least one of the two circles, then this circle is normal. The ordering is trivial resorting to the radii of circles and upper/lower status of θ-monotone circular arcs.

-If the circular arcs intersect transversally, we use the sign of ∆ in Eq. ( 8).

-If the circular arcs are tangent, for k = i or j, we define z k to be either the z-coordinate of the θextremal points of C 0k if it is a normal circle, or the z-coordinate of the pole C 0k goes through if it is a polar circle. We have two sub-cases:

(i)If one of the two circles is threaded, say C 0i , we conclude from the sign of c 0iz -c 0j z if C 0j is threaded too, and from the sign of p z -z j otherwise.

(ii)If both circles are not threaded, we conclude from the radii of the circles and the sign of p z -z k , k = i, j.

Note that if p z -z k is negative (respectively positive), the circular arc of C 0k involved is a lower (respectively upper) one.

Ordering θ-monotone Circular Arcs. We consider two problems: ordering two circular arcs intersecting M θ included in a rational plane, and positioning a point p along a meridian with respect to a circular arc intersected by any given meridian. These two operations are provided by the functor

SK::CompareZAtTheta 3.
For the first predicate, meridian M θ is included in a rational plane. We explicitly construct and compare the z-coordinates of the intersection points between the meridian and the circular arcs. For the second predicate, observe that a non-great circle decomposes S 0 into two regions of unequal areas-the region of largest area and the center of S 0 are on the same side of the plane containing the circle. To begin with, we look for the position of the point p w.r.t. the supporting plane of the circle of the circular arc. This information is sufficient to conclude if p lies on the plane or if the corresponding circle is threaded or polar. When the circular arc lies on a normal circle: If p lies inside the cap of least area defined by the circle on S 0 then p is below the upper circular arc and above the lower one; otherwise the position of p relatively to the circular arc is given by the sign of the difference of p z and the z-coordinate of the θ-extremal points of the circle. Misc. For testing equality of two geometric objects, we compare either their parameters or their algebraic representations. Overlapping tests are performed using underlying geometric objects and inclusion of endpoints. Intersection tests and computations are based on classical mathematical inequalities. Predicate SK::HasOnBoundedSide obviously relies on evaluation of the sign of sphere equation at the point considered. The inclusion tests performed by predicate SK::HasOn 3 are easily triggered using the underlying geometric or algebraic representation. For a point on a circular arc, we additionally need some orientation tests.

Application: Computing Spherical Arrangements

In this section, we present the connexion between the 3D Spherical Kernel and a Bentley-Ottmann like algorithm [START_REF] Bentley | Algorithms for reporting and counting geometric intersections[END_REF] computing the exact arrangement of circles on a sphere [START_REF] Cazals | Computing the exact arrangement of circles on a sphere, with applications in structural biology[END_REF]. In a nutshell, the algorithm takes as input a collection of circles, and returns a decomposition of the sphere into regions whose interiors are connected-the decomposition being stored in an extended half-edge data structure handling holes in faces. In the following, we explain how primitives from the 3D Spherical Kernel are used, and refer the reader to [START_REF] Cazals | Computing the exact arrangement of circles on a sphere, with applications in structural biology[END_REF] for the details on the algorithm.

Algorithm Description. Given a collection of circles on a reference sphere, they are first decomposed into θ-monotone circular arcs. Then, the algorithm consists of sweeping these arcs using a meridian M θ anchored at the poles and moving from 0 to 2π using cylindrical coordinates. This process is tantamount to the classical sweep-line process in the plane, as each θ-monotone circular arc is intersected in at most one point by the sweep meridian. An event corresponds to an intersection point (either a transverse intersection or a tangency point) between two circular arcs, or to a θ-extremal point of a circle. Degeneracies occur when several events are associated to the same point of the reference sphere. To handle them, events are gathered into a data structure called the event site. The vertical ordering V stores the order of θ-monotone circular arcs along the meridian during the sweep. The event queue E stores event sites, which are created upon insertion of intersection and θ-extremal events. Connexion between the Algorithm and the 3D Spherical Kernel. The algorithm features three main constructions. The first two correspond to the construction of θ-extremal points and the creation of θ-monotone circular arcs from input circles using SK::ThetaExtremePoint 3 and SK::MakeThetaMonotone 3 respectively. The third one is concerned with the construction of intersection points and is achieved using

SK::Intersect 3.
Predicates are involved in the manipulation of V and E. Let us first examine the vertical ordering V: (i)The initialization of V requires comparing the position of two circular arcs along meridian M θ with θ = 0, using predicate SK::CompareZAtTheta 3. For intersections occurring on M 0 , one further needs to compare the circular arcs to the right of this point using predicate SK:: CompareZToRight 3. (ii)To update V, the only predicate involved is that required to insert the circular arcs of a normal circle starting to be swept at a given θ-extremal point p. To do so, we locate p along the meridian amongst circular arcs present in V, the predicate involved being SK::CompareZAtTheta 3 with one point and one circular arc. The case of several circles having p as θ-extremal point is easily handled once the positions of the upper and lower arcs of the circle of greatest radius have been determined-this a pencil of circles tangent at p. (iii) To maintain a linear size event queue [START_REF] Brown | Comments on "Algorithms for reporting and counting geometric intersections[END_REF], within an event site, only intersection events corresponding to a pair of circular arcs adjacent along V are stored. This also prevents re-inserting a given event-which is harmful since such a re-insertion causes arithmetic filter failures while seeking this event in the queue. Finally, concatenating appropriately blocks of arcs involved within an event site allows one to maintain V without any numerical operation [START_REF] Cazals | Computing the exact arrangement of circles on a sphere, with applications in structural biology[END_REF].

Let us now analyze the event queue E. Its initialization requires all θ-extremal points. Given all but threaded circles using access function of circles type of circle on reference sphere 3(), such points are constructed using functor SK::ThetaExtremePoint 3. The detection of new intersection points from new adjacencies along V uses predicate SK::DoIntersect 3. All intersection and θ-extremal points are sorted using SK::CompareTheta 3 and SK::CompareZ 3.

The algorithm is summarized in Table 7, the primitives pulled out from the 3D Spherical Kernel being typeset in typewriter font. From this algorithm, the construction of the half-edge data structure storing the arrangement uses two Union-Find processes, as explained in [START_REF] Cazals | Computing the exact arrangement of circles on a sphere, with applications in structural biology[END_REF].

Experiments on Spherical Arrangements

In this section we present two types of results for the algorithm of section 5: On the one hand, we investigate the practical complexity of the algorithm, and on the other hand, we compare to previous work using molecular models.

Variants. The arrangement algorithm being implemented as a generic C++ class, we investigate the following three variants:

-The double variant instantiates the algorithm using a plain double number type as rational type. The main interest is to estimate the overhead imposed by the certification of the results. No guarantee is provided either on the termination or on the correctness of the arrangement.

-The exact variant instantiates the algorithm using CGAL::Lazy Exact NT<Gmpq> as rational type. This is a filtered version of Gmpq [START_REF] Gmp | GNU Multiple Precision Arithmetic Library[END_REF] using the double interval type CGAL::Interval nt. It uses the Gmpq exact rational number type when intervals are not sufficient to certify the answer to a predicate. This number type stores for each number an approximated value and the exact value encoded as a DAG of the arithmetic operations needed to construct it-the value is computed upon request.

-The exact filtered variant instantiates the algorithm twice, respectively with CGAL::Interval nt and CGAL::Lazy Exact NT<Gmpq>, these instantiations being articulated as follows. We first launch the instantiation based upon CGAL::Interval nt. If the certification of a predicate fails, we restart the whole arrangement calculation using the instantiation resorting to CGAL::Lazy Exact NT<Gmpq>. Notice that while the standard filtering strategy consists of re-computing a quantity at the predicate level, we do the same but for a whole arrangement. Notice also that this strategy makes sense if the first calculation fails with low probability. Running Times of an Arrangement of Circles. As shown in [START_REF] Cazals | Computing the exact arrangement of circles on a sphere, with applications in structural biology[END_REF], the Bentley-Ottmann algorithm on a sphere has complexity c(n + k + l) log n, with n the number of circles, k the number of intersection points, and l the number of faces bounded by exactly two arcs (l = O(k) in non-degenerate cases). To scale the running time with the theoretical complexity, we wish to measure to constant c of the algorithm.

To investigate this question, tests were run on collections of random circles on the sphere S 0 centered at the origin, each circle being defined by its center picked uniformly at random within the ball bounded by S 0 . Notice this strategy does not yield great circles. The number of circles chosen varies from 10 to 400 by step size of 10. Computations were run on a bi-proc Pentium IV(R) 3.06Ghz with 2.5GB of RAM.

We performed a curve fitting on running times. As illustrated on Fig. 4, the theoretical curve is in excellent agreement with the experimental curves, for the double and exact variants. Setting the constant c to match the experimental time obtained for n = 400 and the measured k and l, we observe that the constant factor is 2.1 × 10 -5 for the double variant and 10.2 × 10 -5 for the exact variant.

Running Times on Molecular Models. Given a molecular model, we wish to compute one arrangement for each atom-on its surface. This arrangement features the circles defined by the intersections between this atom and its neighbors. These neighbors are retrieved using a regular grid partitioning the axis aligned bounding box of the model into cells. The length of each edge in the grid is the diameter of the largest sphere of the model. Each atomic sphere is associated a cell in the grid and all the spheres intersecting it are either in the same cell or in adjacent ones.

The only running times we are aware of to compute an arrangement of circles on a sphere are those obtained with the algorithm based on explicit controlled perturbations of spheres [START_REF] Eyal | Dynamic maintenance of molecular surfaces under conformational changes[END_REF]. Because the perturbation used is global, the running times reported correspond to the cumulative cost of the arrangements on all spheres. Practically, four protein models from the Protein Data Bank [START_REF]Worldwide Protein Data Bank[END_REF] were used-PDB codes: 1bzm, 1jky, 7at1, 1l7x. Table 8 features our results on these models, using a Pentium III(R) 1Ghz with 1GB of RAM in our case, as opposed to a bi-proc Pentium III(R) 1GHz with 2GB of RAM for [START_REF] Eyal | Dynamic maintenance of molecular surfaces under conformational changes[END_REF]. On these examples, our code is about 65% faster using the double variant and 20% slower using the exact filtered variant.

Consider Tables 8 and9. The ratio between the double variant and the exact filtered variant is of 2. But calculations in double fail on (nearly-)degenerate inputs. One such highly degenerate example is displayed on Fig. 5. Failures also happen for molecular models, since 6 atoms had to be removed to get the double row in Table 9. Input file #atoms From [START_REF] Eyal | Dynamic maintenance of molecular surfaces under conformational changes[END_REF]

Conclusion

High quality geometric code relies on four virtues: robustness, efficiency, modularity, re-usability. Although spheres are amongst the most elementary geometric objects, no library of essential primitives was available to deal with them. This paper answers this need, by developing the cgal 3D Spherical Kernel concept, i.e. a concept featuring the basic types and operations required to deal with spheres, planes, circles, circular arcs and points in 3D. A clear distinction is made between the algebraic and the geometric aspects on the one hand, and on the concepts of a kernel and its implementation on the other hand. The concept is accompanied by an implementation. This implementation, together with a generalization of the Bentley-Ottmann algorithm on a sphere, provides the first effective solution to the problem of computing the exact arrangement of circles on a sphere. Applications in structural biology to investigate protein-protein and protein-drugs interfaces are being developed.

Figure 1 :

 1 Figure 1: The four types of circles on a reference sphere. Black dots are the θ-extremal points.

  represent 3D Spherical Kernel or Algebraic Kernel functors, while names on dashed arrows are access functions.

Figure 3 :

 3 Figure 3: An example of an implementation of a geometric construction: computing the intersection of a circle and a plane.

0. 3 .

 3 Classify circles as normal/polar/bipolar/threaded. SK::CircleORS 3::type of circle on reference sphere 3() Compute θ-extremal points and decompose circles into θ-monotone arcs. SK::ThetaExtremePoint 3, SK::MakeThetaMonotone 3 1. Initialize V: Fill V with circular arcs intersected by the meridian at θ = 0. SK::CompareZAtTheta 3, SK::CompareZToRight 3 2. Initialize E: (a) Look for intersections between circular arcs adjacent in V at θ = 0, SK::DoIntersect 3 and insert the corresponding intersection points into E. SK::Intersect 3, SK::CompareTheta 3, SK::CompareZ 3 (b) For all but threaded circles, insert θ-extremal points into E. SK::CompareTheta 3, SK::CompareZ 3 While E is not empty do (a) Remove from V the circular arcs of the event(s) ending. (b) Insert into V the circular arcs of the event(s) starting.SK::CompareZAtTheta 3 (c) Swap in V the circular arcs intersecting transversally at the event. (d) Insert into E the intersection detected from the new adjacencies along V.SK::DoIntersect 3, SK::Intersect 3, SK::CompareThetaZ 3

Figure 4 :

 4 Figure 4: Observed vs theoretical complexities for random arrangements.

Figure 5 :

 5 Figure5: Left: degenerate arrangement of 47 circles on S 0 . Right: zoom about the point (r 0 , 0, 0): six different directions of tangency for thirteen circles. This point is a θ-extremal point for two circles. Ten circles (two among them being polar) intersect at that point and a bipolar circle passes through it. This arrangement features 674 vertices, 1384 edges, 712 faces and no holes. Euler characteristic reads as 674 -1384 + 712 = 2.

Table 1 :

 1 Sign of δ 3 from Eq. (

	δ 31 δ 32	δ 3	sphere configuration
	+	+	+	S 0 and S i intersect along a non-degenerate circle
	+ --0	-+ -+	--+ 0	S i inside S 0 or S 0 inside S i , no intersection S i outside S 0 , no intersection impossible S i and S 0 tangent, no inclusion
	0 +	-0	0 0	impossible S i and S 0 tangent, S i inside S 0 or S 0 inside S i
	-	0	0	impossible

Table 3 :

 3 Types of the 3D Spherical Kernel. The three groups of the table read as follows: inherited from the linear kernel; general objects; objects on a reference sphere. Notice that objects from the third group inherit the functionalities of objects of the second group.

	Predicates	Returns	Arg.	Meaning
	SK::CompareX 3	{<,=,>}	(p i , p j )	p i .x() vs. p j .x()
	SK::CompareY 3	{<,=,>}	(p i , p j )	p i .y() vs. p j .y()
	SK::CompareZ 3	{<,=,>}	(p i , p j )	p i .z() vs. p j .z()
	SK::CompareXY 3	{<,=,>}	(p i , p j )	Lexicographic order
	SK::CompareXYZ 3	{<,=,>}	(p i , p j )	Lexicographic order
	SK::Equal 3	bool	(w i , w j )	Geometric equality
	SK::HasOn 3 SK::HasOnBoundedSide SK::CompareZAtXY 3	bool bool {<,=,>}	(w i , w j ) (p, S) (p, P )	w j ⊆ w i p ∈ ball bounded by S p vs. P (a, b, d = 0)
	SK::CompareYAtXZ 3	{<,=,>}	(p, P )	p vs. P (a, c, d = 0)
	SK::CompareXAtYZ 3	{<,=,>}	(p, P )	p vs. P (b, c, d = 0)
	SK::DoOverlap 3	bool	(A i , A j )	Overlap test
	SK::DoIntersect 3	bool	(w i , w j )	Intersection test
	SK::CompareTheta 3	{<,=,>}	(p r i , p r j )	Compare θ-coordinates
			(θ i , θ j )	Compare θ-values
	SK::CompareThetaZ 3	{<,=,>}	(p r i , p r j )	Lexicographic order
	SK::CompareZAtTheta 3	{<,=,>}

Table 7 :

 7 The Bentley-Ottmann algorithm for circles on a sphere, together with the primitives of the 3D Spherical Kernel involved.

Table 8 :

 8 Tests on 4 macro-molecular structures. Total time (in seconds) for computing the Van der Waals surface, including the perturbation, taken from[START_REF] Eyal | Dynamic maintenance of molecular surfaces under conformational changes[END_REF] vs. total time of computing the arrangement of circles on each atomic sphere for the three variants.

				Double Exact Exact filtered
	1bzm	2049	8.31	4.98	45.51	9.47
	1jky	5734	26.94	15.75	149.44	29.73
	7at1	7169	28.40	17.00	162.70	32.89
	1l7x	12912	54.50	33.00	311.17	64.32
		NT	neighbor	argt	area total
		double	0.11s	2.14s 0.71s 2.96s
		exact	0.38s	21.47s 7.57s 29.42s
		exact filtered	0.23s	4.01s 1.67s 5.91s

Table 9 :

 9 Comparing Number Types (NT) for complex 1acb (2433 atoms): Run-times to report the neighbors, compute the arrangements, and compute the surface areas of the faces on each atomic sphere.

Acknowledgments

The authors wish to thank Sylvain Pion for helpful discussions.

* This work is partially supported by the IST Programme of the 6th Framework Programme of the EU as a STREP (FET Open Scheme) Project under Contract No IST-006413 -ACS (Algorithms for Complex Shapes with certified topology and numerics)http://acs.cs.rug.nl/

Root of a system of 3 polynomials of the previous types

The AK::RootForSpheres 2 3 type is used to represent the coordinates of a point of the 3D Spherical Kernel. Types. As already mentioned in [START_REF] Emiris | Towards an open curved kernel[END_REF], the Algebraic Kernel must provide basic types (polynomials and roots of systems, see Table 6), and basic functionalities on them. In the sequel, we focus on the most important requirements only.

Main Algebraic Predicates and Constructions. Several elementary functors reflecting the 3D Spherical Kernel user interface are provided for predicates, for instance to compare the coordinates of two AK::RootForSpheres 2 3. Moreover, the functor AK::SignAt allows one to evaluate the sign of a polynomial of type AK::PolynomialForSpheres 2 3 or AK::Polynomial 1 3 at a AK::RootForSpheres 2 3.

The main construction is provided by the functor AK::Solve, which solves a zero-dimensional polynomial system featuring polynomials whose types are listed in Table 6. The return value is an iterator over the solution set (possibly empty), each solution being given as a AK::RootForSpheres 2 3.

Implementation Details

In this section, we bridge the gap between the geometric calculations presented in Section 2 and the functionalities of the kernel: We describe the way concepts developed in the previous section are implemented. Predicates SK::CompareTheta 3 and SK::CompareThetaZ 3 compare two θ-coordinates using a partition of the interval (0, 2π] into eight open intervals of length π 4 and eight value kπ 4 with k = 0 . . . 7. The comparison of θ-coordinates of points falling in different intervals is trivial, while that of points falling in the same interval requires comparing tan θ or cot θ. This boils down to comparing two algebraic number of degree two-each constructed as a quotient (Observation 2) of two degree two algebraic numbers in the same extension using the x-and y-coordinates of each point. This strategy is valid even for events at a pole, using (y i , -x i , 0) and (-y i , x i , 0) as fictitious points whose θ-coordinates match the θ-extremal values of a polar or bipolar circle C 0i . Indeed, one can see that a plane containing the poles and these two points is tangent to or contains C 0i . Identifying Intersection Points. The intersection points of two circles and the θ-extremal values of a circle come into pairs, so that finding the element of the pair with smallest θ-coordinate is an important primitive. Using the afore-mentioned partition of (0, 2π], this is trivial if the two points fall within different intervals. If not, let p and q be two such points. In such a setting, identifying the smallest θ-coordinate, is equivalent to computing Sign(p y q x -p x q y ). By Observation 1 and Corollary 1, there exist four rational numbers a, b, c, d and a rational polynomial of degree two P such that: p x = aR 1 + b p y = cR 1 + d, with P (R 1 ) = 0 q x = aR 2 + b q y = cR 2 + d, with P (R 2 ) = 0 [START_REF] De Castro | Design of the CGAL spherical kernel and application to arrangements of circles on a sphere[END_REF] Therefore, Sign(p y q x -p x q y ) = Sign((bc -ad)(R 1 -R 2 )). Since R 1 and R 2 are roots of the same polynomial, upon creation of p and q we only have to evaluate Sign(bc -ad).

Operations Involving