# Sign of integer determinant 4 by 4

I'm in the context of this publication: http://www.gilbertbernstein.com/resources/booleans2009.pdf

I applied quantization to my point coordinates: All coordinates are integer lying in [0, 2 power 20]. The C++ type representing a coordinate is int64_t.

A plane p is represented in integer homogeneous coordinates (a, b, c, d):

p = { (x, y, z) | ax + by + cz + d = 0 }


A point X can be represented as 3 planes (p, q, r) (provided a validity condition, cf the publication).

I have a valid point X and a plane s. I want to know if X lies on the positive side, negative side or on s. According to the publication, I need to compute the sign of these determinants:

| p_a, p_b, p_c, p_d |
| q_a, q_b, p_c, q_d |
| r_a, r_b, p_c, r_d |
| s_a, s_b, p_c, s_d |


and:

| p_a, p_b, p_c|
| q_a, q_b, p_c|
| r_a, r_b, p_c|


The quantization bound implies the following bounds:

|p_a|, |p_b| and |p_c| < 2 pow 40


and

|p_d| < 2 pow 60


How to compute exactly and quickly the sign of that determinant? Using BigInt or anything based on the Chinese Remainder theorem is too slow (I haven't tried the GMP library though)... What libraries that could solve this problem?

Clue:

There are several options for deriving predicates with arbitrary precision.

If you use integers, then you can use Haramard inequality (see standard algebra textbooks, for instance Greub 1975, Linear Algebra, Springer Verlag). This inequality says that the absolute value of a determinant is not greater than the product of the lengths of the column vectors. You can use it to bound the number of bits of the determinant in function of the number of bits of the coefficients. For the 4x4 determinant of the in_circle predicate of a 2D Voronoi diagram, if coordinates are between $-K$ and $K$, the absolute value of the determinant is smaller than $32K^4$.If you use long integers (64 bits), for a 2D Voronoi diagram, $K=23170$ ensures that you got sufficient precision. If you are in 3D, then 64 bits may not suffice, and you will need to use a bigint package (GNU GMP as you suggested). It is probably possible to first do the computation with standard 64 bit integers and detect dynamically whether an overflow occured, then relaunch the computation with the bigint package only when needed.

If you use floating point numbers, see my answer to this question, with reference to my survey article and software package. The idea is to compute the determinant with floating point numbers as well as a bound indicating whether the sign is accurate, then relaunch the computation with (more expensive) arbitrary precision if this was not the case. My software package 2 can be used to derive predicates involved in intersections, as explained in my article 3. This uses Shewchuk's expansions (mentioned in the question) combined with an arithmetic filter 4. This strategy is easier to implement than the adaptive evaluation used in Shewchuk's original package, and it is reasonably easy to derive the new predicates that you will need for your intersections.

more thoughts / references: I recommend the following article by a friend of mine (QuickCSG) on the very same topic (robust evaluation of boolean operations on meshes). It recasts the intersection problem as an algorithm that only depends on a very small number of predicates, for which valid simulation of simplicity can be reasonably easily defined.

• Thanks a lot Bruno, I'm extremely pleased by such an informative and deep answer! the floating point numbers solution seems the most promising, I'll keep you updated about my results. – Brainless Sep 11 '17 at 20:58

Essentially, we are trying to compute the determinant of integer matrices. See if these helps:

Exact Determinant of Integer Matrices, Takeshi Ogita http://www.eng.nus.edu.sg/civil/REC2010/documents/papers/038.pdf

Computing the sign or the value of the determinant of an integer matrix: A complexity survey, Erich Kaltofena, Gilles Villard http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.21.9676&rep=rep1&type=pdf

Also, if we think in more geometric terms, we could always assemble a numerically more stable, easy to compute way to determine this. Let $\mathbf{n}=(a,b,c)$ denote the normal of the plane and $\mathbf{p} \in \mathbb{R}^3$ the point. Let's take a random point on this plane $\mathbf{p}_0=[0,0,-c/d]^T$ and construct a direction vector $\mathbf{d}=\mathbf{p}-\mathbf{p}_0$. The point lies on the positive side of the plane iff: $$\mathbf{n}\cdot\mathbf{d} > 0 \rightarrow \text{positive side}\\ \mathbf{n}\cdot\mathbf{d} < 0 \rightarrow \text{negative side}$$

doing these operations might be fairly easy also for integer valued/quantized input, where BigInt can actually be used.

• Thanks for your answer Tolga, unfortunately since we're in discrete geom (coordinates are integers), the coordinates of a point represented as the intersections of 3 planes might not have integer coordinates... And the simple intersection test you wrote may not be accurate. Furthermore, computing the determinant, even with the papers you mentioned, will be too slow I think... – Brainless Sep 11 '17 at 21:01
• But any proposition is still welcome! – Brainless Sep 11 '17 at 21:01