Certified geometric predicates

For each row $i$, write its entries as $M_{ij}=m_{ij}{2}^{e_{ij}}$ using [F1], and choose $k_i=\max (0,-{\min }_j{e}_{ij})$. The row is finite, so this integer exists; then every $N_{ij}=2^{k_i}{M}_{ij}$ is an integer.

In every Leibniz term of $\det N$, the product contains exactly one entry from each row, so it equals $({\prod }_i{2}^{k_i})$ times the corresponding term of $\det M$. Summing the finite terms using [F2] gives $\det N=({\prod }_i{2}^{k_i})\det M$.

Every factor $2^{k_i}$ is positive. Therefore $\det N$ and $\det M$ have the same sign and vanish simultaneously.

For an integer $n\ge 0$ with $nu<1$, the binomial theorem and $\left(\genfrac{}{}{0px}{}{n}{k}\right)\le n^k$ give $(1+u{)}^n\le {\sum }_{k\ge 0}(nu{)}^k=1\mathrm{/}(1-nu)$. Therefore $(1+u{)}^n-1\le {\gamma }_n$.

Put $a=s-1$ and let $t_{\pi }$ be one exact signed Leibniz product. Its $a$ rounded multiplications introduce factors $(1+{\delta }_k)$ by [F1] and [A1]. The identity ${\prod }_k(1+{\delta }_k)-1={\sum }_k{\delta }_k{\prod }_{j<k}(1+{\delta }_j)$, the triangle inequality [F2], and $\mathrm{\mid }1+{\delta }_j\mathrm{\mid }\le 1+u$ show $\mathrm{\mid }{\widehat{t}}_{\pi }-t_{\pi }\mathrm{\mid }\le ((1+u{)}^a-1)\mathrm{\mid }{t}_{\pi }\mathrm{\mid }\le {\gamma }_a\mathrm{\mid }{t}_{\pi }\mathrm{\mid }$ by [step 1.1]. The same product bound gives $\mathrm{\mid }{\widehat{t}}_{\pi }\mathrm{\mid }\le (1+{\gamma }_a)\mathrm{\mid }{t}_{\pi }\mathrm{\mid }$; an exact zero term remains zero.

There are $s!$ terms and hence $b=s!-1$ rounded additions after the first term is placed in the accumulator. If the successive addition errors are ${\eta }_1,\dots ,{\eta }_b$, induction in the accumulator recurrence writes $\widehat{D}={\sum }_{j=1}^{s!}{\widehat{t}}_j{P}_j$, where each $P_j$ is a product of at most $b$ factors $(1+{\eta }_k)$. The telescoping identity in [step 2.1] gives $\mathrm{\mid }{P}_j-1\mathrm{\mid }\le {\gamma }_b$. Thus $\mathrm{\mid }\widehat{D}-{\sum }_j{\widehat{t}}_j\mathrm{\mid }=\mathrm{\mid }{\sum }_j{\widehat{t}}_j(P_j-1)\mathrm{\mid }$ is at most ${\gamma }_b{\sum }_j\mathrm{\mid }{\widehat{t}}_j\mathrm{\mid }$ by [F2].

Apply [F2], [step 2.1], and [step 3.1] to the exact signed sum. With $S={\sum }_{\pi }\mathrm{\mid }{t}_{\pi }\mathrm{\mid }$ this gives $\mathrm{\mid }\widehat{D}-\det M\mathrm{\mid }\le ({\gamma }_a+{\gamma }_b+{\gamma }_a{\gamma }_b)S$.

Since $q=a+b$ and $qu<1$, $(1-au)(1-bu)=1-qu+ab{u}^2\ge 1-qu>0$. Therefore $(1+{\gamma }_a)(1+{\gamma }_b)-1\le {\gamma }_q$, so the coefficient in [step 4.1] is at most ${\gamma }_q$. This proves the stated bound.

If $d=2$, write $A=\left(\begin{array}{cc}a& b_0\\ c& d_0\end{array}\right)$ and $\mathrm{\Delta }=a{d}_0-b_{0}c=\det A$ by [F1]. For $B=\left(\begin{array}{cc}{d}_0& -b_0\\ -c& a\end{array}\right)$, entrywise multiplication gives $AB=\left(\begin{array}{cc}a{d}_0-b_{0}c& -a{b}_0+a{b}_0\\ c{d}_0-d_{0}c& -c{b}_0+d_{0}a\end{array}\right)=\mathrm{\Delta }{I}_2$ and $BA=\left(\begin{array}{cc}{d}_{0}a-b_{0}c& d_0{b}_0-b_0{d}_0\\ -ca+ac& -c{b}_0+a{d}_0\end{array}\right)=\mathrm{\Delta }{I}_2$. Since $\mathrm{\Delta }\mathrm{\ne }0$, [F2] gives ${\mathrm{\Delta }}^{-1}$, and ${\mathrm{\Delta }}^{-1}B$ is a two-sided inverse of $A$.

If $d=3$, write $A=\left(\begin{array}{ccc}a& b_0& c\\ d_0& e& f\\ g& h& i\end{array}\right)$ and $\mathrm{\Delta }=aei+b_{0}fg+c{d}_{0}h-ceg-b_0{d}_{0}i-afh=\det A$ by [F1]. Set $B=\left(\begin{array}{ccc}ei-fh& ch-b_{0}i& b_{0}f-ce\\ fg-d_{0}i& ai-cg& c{d}_0-af\\ d_{0}h-eg& b_{0}g-ah& ae-b_0{d}_0\end{array}\right)$. Direct multiplication gives diagonal entries of $AB$ equal respectively to $aei-afh+b_{0}fg-b_0{d}_{0}i+c{d}_{0}h-ceg$, $c{d}_{0}h-b_0{d}_{0}i+aei-ceg+b_{0}fg-afh$, and $b_{0}fg-ceg+c{d}_{0}h-afh+aei-b_0{d}_{0}i$, each $\mathrm{\Delta }$. Its off-diagonal entries, in row order, are $ach-a{b}_{0}i+a{b}_{0}i-b_{0}cg+b_{0}cg-ach=0$, $a{b}_{0}f-ace+c{b}_0{d}_0-a{b}_{0}f+ace-c{b}_0{d}_0=0$, $d_{0}ei-d_{0}fh+efg-d_{0}ei+d_{0}fh-efg=0$, $d_0{b}_{0}f-d_{0}ce+ec{d}_0-eaf+fae-f{b}_0{d}_0=0$, $gei-gfh+hfg-h{d}_{0}i+i{d}_{0}h-ieg=0$, and $gch-g{b}_{0}i+hai-hcg+i{b}_{0}g-iah=0$. Thus the displayed cofactor formula, denoted $B(C)$ when applied to any $3\times 3$ matrix $C$, satisfies $CB(C)=\det (C)I_3$ by the same entrywise calculation. Substitution of $A^T$ into that formula gives $B(A^T)=\left(\begin{array}{ccc}ei-hf& gf-d_{0}i& d_{0}h-ge\\ hc-b_{0}i& ai-gc& g{b}_0-ah\\ b_{0}f-ec& d_{0}c-af& ae-d_0{b}_0\end{array}\right)=B^T,$ where the equality uses commutativity of the real field. In the Leibniz sum, reindexing $\pi$ by $\sigma ={\pi }^{-1}$ gives $\det (A^T)=\det A$: the product ${\prod }_i{A}_{\pi (i),i}$ becomes ${\prod }_i{A}_{i,\sigma (i)}$, and $\mathrm{sgn}(\pi )=\mathrm{sgn}(\sigma )$ because the inversion map $(i,j)\mapsto (\pi (j),\pi (i))$ bijects inversions of $\pi$ with inversions of ${\pi }^{-1}$. Therefore applying the displayed cofactor identity to $A^T$ yields $A^T{B}^T=\det (A^T)I_3=\mathrm{\Delta }{I}_3$. Transposing this equality reverses the factors by the finite matrix-product formula, so $BA=\mathrm{\Delta }{I}_3$. Since $\mathrm{\Delta }\mathrm{\ne }0$, ${\mathrm{\Delta }}^{-1}B$ is a two-sided inverse of $A$ by [F2].

In either dimension, let $C$ be the two-sided inverse obtained in [step 1.1] or [step 1.2]. Then $x=Cb$ satisfies $Ax=(AC)b=I_{d}b=b$. If $Ax=Ay$, multiplying by $C$ gives $x=(CA)x=(CA)y=y$. Hence the solution exists and is unique.

Let $A$ have row $i$ equal to $p_i-p_0$ for $i=1,\dots ,d$. Subtracting the first row from every later row of the matrix defining $O_d$ preserves its determinant: multilinearity in each row follows from the finite Leibniz sum, and every added determinant has two equal rows, whose permutation terms cancel in pairs. Expanding along the final column gives $O_d=(-1{)}^d\det A$, so $\det A\mathrm{\ne }0$.

Put $b_i=({\sum }_j{p}_{i,j}^2-{\sum }_j{p}_{0,j}^2)\mathrm{/}2$ for $i=1,\dots ,d$. Since $1>0$ and the order is compatible with addition, $2=1+1>0$ and is nonzero. By [F3] and [step 1.1], there is a unique $c\in {\mathbb{R}}^d$ with $Ac=b$.

Expanding $Ac=b$ gives $2{\sum }_j{c}_j(p_{i,j}-p_{0,j})={\sum }_j{p}_{i,j}^2-{\sum }_j{p}_{0,j}^2$. Rearranging shows ${\sum }_j(p_{i,j}-c_j{)}^2={\sum }_j(p_{0,j}-c_j{)}^2$ for every $i=1,\dots ,d$. Let $s={\sum }_j(p_{0,j}-c_j{)}^2\ge 0$; by [F4], the unique $r\ge 0$ with $r^2=s$ exists. Hence $S(c,r)$ passes through all base points.

If $(c^{\mathrm{\prime }},r^{\mathrm{\prime }})$ is any other center and nonnegative radius through the base points, equal squared distances give $A{c}^{\mathrm{\prime }}=b$. Uniqueness in [F3] forces $c^{\mathrm{\prime }}=c$. Then $r^{\mathrm{\prime }2}=s=r^2$, and uniqueness of the nonnegative square root [F4] gives $r^{\mathrm{\prime }}=r$. All objects are obtained from finite data and unique solutions; no choice axiom is used.

For each admissible choice of coordinates inside the input intervals, the determinant expression in [F1] is an exact evaluation of the finite DAG used by the algorithm. If all DAG operations are finite, [F2] places that exact determinant in the computed interval $J$.

If the procedure returns positive, [F3] gives $J=[a,b]$ with $a>0$; every admissible determinant in $J$ is positive. If it returns negative, $b<0$ and every admissible determinant is negative. When $0\in J$ or a nonfinite status occurs it returns indeterminate and makes no uniform-sign claim.

Let $A$ be the $d\times d$ matrix with rows $p_i-p_0$ for $i=1,\dots ,d$. Subtracting row $p_0$ from the other rows of the matrix defining $O_d$ leaves its determinant unchanged: the Leibniz formula is linear in each row, and the added determinant has two equal rows, whose terms cancel in pairs under the transposition of those rows. Expansion along the final column then gives $\det A=(-1{)}^d{O}_d$.

Expanding the equal sphere radii in [A1], which exist uniquely by [F3], gives ${\sum }_j{p}_{i,j}^2-{\sum }_j{p}_{0,j}^2=2{\sum }_j{c}_j(p_{i,j}-p_{0,j})$ for $i=1,\dots ,d$. For $q$, put $R={\sum }_j(q_j-c_j{)}^2-r^2$; the same expansion gives ${\sum }_j{q}_j^2-{\sum }_j{p}_{0,j}^2=2{\sum }_j{c}_j(q_j-p_{0,j})+R$.

Subtract row $p_0$ from the other rows of the matrix defining $L_d$, then expand along its final column. The cofactor sign is $(-1{)}^{d+3}$, and the remaining determinant has rows $(p_i-p_0,{\sum }_j{p}_{i,j}^2-{\sum }_j{p}_{0,j}^2)$, including the row for $q$. Subtracting $2c_j$ times coordinate column $j$ from the last column makes its first $d$ entries zero by [step 1.2] and its final entry $R$. The resulting determinant is $\det AR$. These row and column operations preserve the Leibniz determinant: multilinearity in columns follows from the same permutation formula, and the added determinant has equal columns whose terms cancel in pairs under their transposition. Thus $L_d=(-1{)}^{d+3}\det AR=-O_{d}R$.

Multiplication by $O_d$ gives $O_d{L}_d=-O_d^{2}R$. Since $O_d\mathrm{\ne }0$, $O_d^2>0$ in the ordered field [F2]; therefore the product is positive exactly when $R<0$, zero exactly when $R=0$, and negative exactly when $R>0$.

The filter's absolute-term DAG evaluates $S$ from [F2]. Its returned upper endpoint $U$ satisfies $S\le U$ by [F3]. The definition of $E$ and [F4] give $E\ge {\gamma }_{q}U$, hence $E\ge {\gamma }_{q}S$.

Combining [F2] with [step 1.1] yields $\mathrm{\mid }\widehat{D}-\det M\mathrm{\mid }\le E$, so $\det M\in [\widehat{D}-E,\widehat{D}+E]$.

If the algorithm returns positive, [F1] says $\widehat{D}-E>0$; by [step 2.1], $\det M\ge \widehat{D}-E>0$. If it returns negative, then $\widehat{D}+E<0$ and [step 2.1] gives $\det M\le \widehat{D}+E<0$. In all other cases it returns no sign.

Each finitely encoded dyadic entry can be put in the form in [F1] by integer shifts and cancellation of powers of two; zero uses exponent zero. Since $s\ge 1$, every row has a nonempty finite list of exponents and its maximum $k_i$ exists and is computed by finitely many comparisons. Then $N_{ij}=a_{ij}{2}^{k_i-b_{ij}}$ is an exactly computable integer.

In each permutation product of $N$, the row factors multiply to the same positive integer $C={\prod }_i{2}^{k_i}$. Thus [F2] gives $\det N=C\det M$, the sign-preserving identity underlying [F3]. The integer sum in [F1] is exactly $\det N$, so the returned sign equals $\mathrm{sgn}(\det M)$, also when this value is zero.

There are exactly $s!$ permutations, which can be enumerated in lexicographic order from the finite set of tuples of length $s$ with entries in $\{0,\dots ,s-1\}$ by retaining those with distinct entries. Each product, sum, integer shift and sign comparison operates on finite integer encodings and terminates. Hence all loops and the algorithm terminate, even for singular matrices and $s=1$; there is no pivot division. All constructions are explicit and finite, so no axiom of choice is used.

Regard one requested determinant as a polynomial in independent coordinate increments $z_{i,j}$. It is not the zero polynomial: for orientation choose the perturbed points $0,e_1,\dots ,e_d$, whose homogeneous determinant is nonzero. For a lifted determinant choose the first $d+1$ perturbed points $0,e_1,\dots ,e_d$ and the last point $(2,\dots ,2)$. The base sphere has center $(1\mathrm{/}2,\dots ,1\mathrm{/}2)$ and squared radius $d\mathrm{/}4$, while the last point has squared distance $9d\mathrm{/}4$ from that center; [F3] therefore gives a nonzero lifted determinant. These coordinate assignments are available because the $z_{i,j}$ are independent variables.

In every orientation determinant each coordinate increment has degree at most one. In a lifted determinant it has degree at most two, because it can occur squared in one lifted entry and a determinant term selects only one entry from that row. Thus every monomial exponent vector has digits $a_k\in \{0,1,2\}$. Under the substitution in [F2], its $\epsilon$ exponent is ${\sum }_k{a}_k{3}^k$. Uniqueness of base-three expansion makes this map injective, so distinct monomials cannot merge or cancel. The nonzero polynomial from [step 1.1] remains nonzero after substitution.

Write any resulting polynomial as $P(\epsilon )=c_m{\epsilon }^m+{\sum }_{j>m}{c}_j{\epsilon }^j$, with $c_m\mathrm{\ne }0$. If $B={\sum }_{j>m}\mathrm{\mid }{c}_j\mathrm{\mid }$ is positive, then for $0<\epsilon <\min (1,\mathrm{\mid }{c}_m\mathrm{\mid }\mathrm{/}(2B))$ the tail after dividing by ${\epsilon }^m$ has absolute value at most $\epsilon B<\mathrm{\mid }{c}_m\mathrm{\mid }\mathrm{/}2$ by [F5] and [F6]. If $B=0$, its sign is already constant. In either case the sign for sufficiently small positive $\epsilon$ is the sign of $c_m$. If the original determinant is nonzero, its constant coefficient is nonzero and is the first coefficient, so its sign is preserved.

A finite point set has only finitely many ordered orientation and lifted predicate tuples. For each resulting polynomial with tail bound $B>0$ in [step 3.1], form its positive threshold $t_P=\min (1,\mathrm{\mid }{c}_m\mathrm{\mid }\mathrm{/}(2B))$. If the finite set of these thresholds is nonempty, let $t$ be its minimum and choose $\epsilon =t\mathrm{/}2$; then [step 3.1] applies to every such polynomial. Every remaining polynomial has $B=0$, so it is a single nonzero monomial and has its first-coefficient sign for every $\epsilon >0$. If there are no thresholds, choose $\epsilon =1$; the same reasoning covers the case of an empty predicate family. Thus one positive value realizes all returned signs simultaneously and no declared determinant is zero. The algorithm returns exactly these first-coefficient signs without selecting or evaluating that $\epsilon$, so the tie-break is consistent and uses no arbitrary epsilon threshold.

The point tuple is finite and every coordinate is dyadic. The entries of the orientation matrix and lifted matrix use only finitely many sums and products of dyadics, so each requested matrix is a finite dyadic matrix under [F1].

For each requested matrix, if the floating filter returns a sign, it is exact by [F2]. If the filter is skipped or returns inconclusive, the algorithm invokes the integer fallback; [F3] computes its finite Leibniz sum exactly, and [F4] transfers its sign to the original dyadic determinant. Thus the raw result is correct, including when it is zero.

A matrix has size at most five. Each filter attempt enumerates finitely many products and interval operations, and a failed precondition or nonfinite status branches directly to the fallback. The fallback performs finite row scalings and exactly $s!$ integer terms, each a finite integer product and sum, so it terminates for every finite encoding.

When a sphere-side result is requested and the base orientation is nonzero, the algorithm returns the sign of the product of the two exact raw signs from [step 1.2]. Its inside/on/outside meaning follows from [F5]. When the base orientation is zero, the algorithm returns degenerate-base and makes no side classification.

The construction and fallback are finite and choose no element from an arbitrary family; no axiom of choice is used.