Rational Points on the Unit Sphere: Approximation Complexity and Practical Constructions
Daniel Bahrdt Martin P. Seybold
cs.CG

	
	
	Improved Analysis<cit.>
	
	
	
	Rational Points on the Unit Sphere: Approximation Complexity and Practical Constructions
     Daniel Bahrdt	Formale Methoden der Informatik, University of Stuttgart,  Germany, {bahrdt, seybold}@fmi.uni-stuttgart.de Martin P.  Seybold[1]
    December 30, 2023
===================================================================================================================================================



  §.§.§ Abstract
 	

		Each non-zero point in ^d identifies a closest point x on the unit sphere ^d-1.
		We are interested in computing an ϵ-approximation y ∈^d for x, that is exactly on ^d-1 and has low bit size. We revise lower bounds on rational approximations and provide explicit, spherical instances.
		
		We prove that floating-point numbers can only provide trivial solutions to the sphere equation in ^2 and ^3. Moreover, we show how to construct a rational point with denominators of at most 10(d-1)/ε^2 for any given ϵ∈(0, 1 8], improving on a previous result.
		The method further benefits from algorithms for simultaneous Diophantine approximation.
		
		Our open-source implementation and experiments demonstrate the practicality of our approach in the context of massive data sets Geo-referenced by latitude and longitude values.
	

	

  
*Keywords
	 Diophantine approximation,
		Rational points,  
		Unit sphere,
		Perturbation, 
		Stable geometric constructions
	
	





§ INTRODUCTION

Many mathematical sciences use trigonometric functions in symbolic coordinate transformations to simplify fundamental equations of physics or mathematical systems.
However, rational numbers are dominating in computer processing as they allow for simple storage as well as fast exact and inexact arithmetics (e.g. GMP<cit.>, IEEE Float, MPFR<cit.>).
Therefore problems on spherical surfaces often require to scale a point vector, as in choosing a point  uniform at random<cit.>, or to evaluate a trigonometric function for a rational angle argument, as in dealing with Geo-referenced data.

A classical theoretical barrier is Niven's theorem<cit.>, 
 which states that the sole rational values of sine for rational multiplies of π are 0, ±1/2 and ± 1.
The well known Chebyshev polynomials have roots at these values, hence give rise to representations for these algebraic numbers.
However, arithmetics in a full algebraic number field might well be too demanding for many applications.
For products of sine and cosine, working with Euler's formula on the complex unit circle and Chebyshev polynomials would suffice though.

This manifests in problems of exact geometrical computations, since standard methodology relies on Cartesian input<cit.>.
Spheres and ellipsoids are common geometric objects and rational solutions to their defining quadratic polynomials are closely related to Diophantine equations of degree 2.
The famous Pythagorean Triples are known to identify the rational points on the circle ^1.
Moreover, the unit sphere has a dense set of rational points and so do ellipsoids with rational half-axes through scaling.
Spherical coordinates are convenient to reference such Cartesians with angle coordinates and geo-referenced data denotes points with rational angles.
Standard approximations of Cartesians do not necessarily fulfill these equations, therefore subsequent algorithmic results can suffer greatly.	
	
This paper focuses on finding rational points exactly on the unit sphere ^d-1={x ∈^d  : ∑_i x_i^2 = 1} with bounded distance to the point x / ‖ x ‖_2 – its closest point on ^d-1.
In this work, x ∈^d can be given by any finite means that allow to compute a rational approximation to it with arbitrary target precision.
Using rational Cartesian approximations for spherical coordinates, as derived from MPFR, is just one example of such a black-box model.
Moreover, we are interested in calculating rational points on ^d with small denominators.





 §.§ Related Work
 
Studies on spherical Delaunay triangulations (SDT), using great-circle segments on the sphere ^2, provide common ways to avoid and deal with the point-on-sphere problem in computational geometry.

The fragile approaches
<cit.>
ignore that the input may not be on ^2 and succeed if the results of all predicate evaluations happen to be correct.
Input point arrangements with close proximity or unfortunate locations bring these algorithms to crash, loop or produce erroneous output.

The quasi-robust approaches <cit.> weaken the objective and calculate a Delaunay tessellation in d-Simplexes.
Lifting to a d+1 convex hull problem is achieved by augmenting a rational coordinate from a quadratic form – The augmented point exactly meets the (elliptic) paraboloid equation.
However, the output only identifies a SDT if all input points are already on the sphere, otherwise the objectives are distinct. Equally unclear is how to address spherical predicates and spherical constructions.

The robust approaches <cit.> use the circle preserving stereographic projection from ^2 to the plane.
The perturbation to input, for which the output is correct, can be very large as the projection does not preserve distances. 
Furthermore, achieving additional predicates and constructions remains unclear.

The stable approaches provide geometric predicates and constructions for points on ^2 by explicitly storing an algebraic number, originating from scaling an ordinary rational approximation to unit length<cit.>.
Algebraic number arithmetics can be avoided for ^2, but exact evaluation relies on specifically tailored predicates <cit.>, leaving the implementation of new constructions and predicates open.



Kleinbock and Merrill provide methods to quantify the density of rational points on ^d <cit.>, that extend to other manifolds as well.
Recently, Schmutz<cit.> provided an divide-&-conquer approach on the sphere equation, using Diophantine approximation by continued fractions, to derive points in ^d ∩^d-1 for a point on the unit sphere ^d-1.
The main theorem bounds the denominators in ε-approximations, under the ‖ ‖_∞ norm, with
( √(32)⌈log_2 d ⌉ / ε  )^2⌈log_2 d ⌉.
Based on this, rational approximations in the orthogonal group O(n,) and in the unitary matrix group U(n,) are found. 
This is of particular interest for sweep-line algorithms:
<cit.> studies finding a rotation matrix with small rationals for a given rational rotation angle of an 2D arrangement.




 §.§ Contribution

The strong lower bound on rational approximations to other rational values does not hold for Geo-referenced data considering Niven's theorem.
We derive explicit constants to Liouville's lower bound, for a concrete Geo-referenced point, that is within a factor 2 of the strong lower bound.
Moreover, we prove that floating-point numbers cannot represent Cartesian coordinates of points that are exactly on ^1 or ^2.

We describe how the use of rotation symmetry and approximations with fixed-point numbers suffice to improve on the main theorem of <cit.>.
We derive rational points exactly on ^d-1 with denominators of at most
10(d-1)/ε^2
for any ε∈(0, 1 8].
Moreover, our method allows for even smaller denominators based on algorithms for simultaneous Diophantine approximations, though a potentially weaker form of approximation would suffice.

The controlled perturbations provided by our method allow exact geometric algorithms on ^d to rely on rational rather than algebraic numbers – E.g. enabling convex hull algorithms to efficiently obtain spherical Delaunay triangulations on ^d and not just Delaunay tessellations.
Moreover, the approach allows for inexact but ε-stable geometric constructions – E.g. intersections of Great Circle segments.

We demonstrate the quality and effectiveness of the method on several, including one whole-world sized, point sets.
We provide open-source implementations for the method and its application in the case of spherical Delaunay triangulations with intersections of constraints.



§ DEFINITIONS AND TOOLS

The 2nd Chebyshev polynomials U_n of degree n are in [X], given their recursive definition:

    U_0(x)   =1      U_1(x)=2x
    
    U_n+1(x)   =2xU_n(x)-U_n-1(x)  .

It is well known <cit.>, that the n roots of U_n are exactly the values  
{ cos(π k / (n+1)) : k=1,… ,n  }.
Hence the polynomials U_n give rise to algebraic representations for cosine values of rational multiplies of π.
This is particularly useful in conjunction with classic results on Diophantine approximations, that are known since 1844<cit.>:

 
	For any algebraic α∈ of degree n ≥ 2,
there is a positive constant c(α)> 0 such that 
	
    | α - p/q| ≥c(α)/q^n

for any p ∈ and q ∈.

Apart from this lower bound on rational approximations, there is another important folklore result on the existence of simultaneous Diophantine approximations.
Such approximations have surprisingly small errors, despite their rather small common denominator.
 
	Let N ∈ and α∈^d with 0 ≤α_i ≤ 1.
	There are integers p ∈^d,  q ∈ with 1 ≤ q ≤ N and 
	
    |α_i - p_iq|≤1/q√(N) .


	(See Appendix <ref> for a proof.)
For d=1, the continued fraction (equivalently the Euclidean) algorithm is famous <cit.> for finding approximations with | α - p/q| ≤ 1/2q^2.
This spurred the field of number theory to study generalizations of the continued fraction algorithm that come close to Dirichlet's upper bound, but avoid brute-force calculations.
Some more recent methods are discussed in Section <ref>.

Our approach uses the Stereographic Projection in ^d.
Let p=(0,…,0,1) ∈^d be the fixed point for the projection τ, mapping all points of a ray from p to the intersection with the hyperplane x_d = 0.

    τ : ^d∖(^d-1×{1})    →^d-1
    
    x    ↦(x_1/1-x_d , … , x_d-1/1-x_d)

The surjective mapping τ is injective as well, when restricted to the domain ^d-1∖{p}.
We further define the mapping σ, which is

    σ : ^d-1   →^d∖{p}
    
    x    ↦(  2x_1/1+S^2 , … , 2x_d-1/1+S^2 , -1+S^2/1+S^2 )

where S^2 = ∑_j=1^d-1x_j^2.
We have  imgσ⊆^d-1, since 

    ‖σ(x) ‖_2^2= (-1+S^2)^2 + ∑_i=1^d-1 (2x_i)^2/(1+S^2)^2 = 1  .

Furthermore, x = τ∘σ( x) for all x∈^d-1, since

    (τ∘σ)_i(x) =
    2x_i/1+S^2/1--1+S^2/1+S^2 = 2x_i/1+S^2+1-S^2 =x_i

holds for all 1 ≤ i < d.
Hence, σ and τ are inverse mappings.
Note that images of rational points remain rational in both mappings, establishing a bijection between rational points in ^d-1 and ^d-1.




 §.§ Lower Bounds and Instances for Geo-referenced Data on ^d
 
It is well known in Diophantine approximation that rational numbers have algebraic degree 1 and are hard (in the following qualitative sense) to approximate with other rational numbers.
The following folklore observation is an analog to Liouville's lower bound.


	For  rational numbers ab ≠ pq, we have
	
    | a/b- p/q| =	| aq-bp/bq| ≥1/bq


If q<b, we have a lower bound of 1/q^2 for rational approximations to ab with denominators up to q.
Pythagorean triples (x,y,z) ∈^3 provide such rational points on ^1, since (x/z)^2 + (y/z)^2 =1.
We have a lower bound of 1/z^2 for approximations with denominators q<z.
See Section <ref> for rational points on ^d with the same denominator property.

The situation might look different when dealing with Geo-re­fer­enced data (rational angle arguments) only.
However, using Chebyshev's polynomials in conjunction with Liouville's lower bound (c.f. Theorem <ref>) allows to derive explicit constants for Diophantine approximations of cos(108^∘).

Given spherical coordinates, the first coordinate of a point on ^d might well have algebraic values of r_i= cos( i5π) for i ∈{1,2,3,4}.

    (r_1,r_2,r_3,r_4)
       =( 1+√(5)/4, -1+√(5)/4, 1-√(5)/4, -1-√(5)/4)
       ≈( +0.8090, +0.3090, -0.3090, -0.8090)

Over [X], the polynomial U_4(x)=16x^4-12x^2+1 has the irreducible factors

    U_4(x) = (4x^2 - 2x + 1)_=:f(x)(4x^2 + 2x - 1)

Since r_1 and r_3 are the roots of the polynomial f, they have algebraic degree n=2.

Using Liouville's lower bound for r_3, we have for all pq ∈

    | r_3 - p/q| ≥min{c_2,  1c_1}/q^n  ,

with constants c_1 and c_2 according to the proof of Liouville's Theo­rem<cit.>.
The constants c_1,c_2>0 exist, since the polynomial division of f with the linear factor (x-r_3) results in the continuous function g(x)=(x-r_1).
For c_2 = 1/2 < √(5)/2, the interval I:=[r_3- c_2, r_3 +c_2] ⊆ is sufficiently small to exclude different roots of f and the inequality

    max_x∈ I|g(x)| = 	max_x∈ I|x - r_1| < c_1

is met with a generous choice of c_1=2.
This leads to an explicit lower bound on the approximation error to r_3 with denominators q of

    | cos(108^∘) - p/q| ≥1/2 · q^2 .




§ RESULTS

Apart from integers, contemporary computing hardware heavily relies on floating point numbers.
These are triplets (s,m,e) with s ∈{0,1}, m ∈{0,…,2^l-1} and e ∈{-2^k-1+1,…, 2^k-1-1}.
The IEEE standard for Float is (l,k)=(23,8) and (52,11) for Double.
The rational number described by such a triplet is 

    val(s,m,e) = (-1)^s ·2^l + m2^l  2^e   							    e > 0  
    
    2^l + m2^l12^|e|    e < 0  
    
      0 + m2^l12^2^k-1-2    e = 0

where the latter case describes `denormalized' numbers.
In each case, the uncanceled rational value has some power of 2 as the denominator.
Since powers of two are the sole divisors of a 2^i, the denominator of the canceled rational has to be a power of two, too.
Hence, rational values representable by floating point numbers are a subset of the following set P and fixed-point binary numbers are a subset of P_i:

    imgval⊆{z/2^i : i ∈, z ∈, z  odd}=    P 
    {z/2^i : z ∈} = P_i ⊆   P   .

	


 §.§ Floating Point Numbers are Insufficient

Fix-point and floating-point arithmetics of modern CPUs work within a subset of rational numbers, in which the denominator is some power of two and the result of each arithmetic operation is `rounded'.

	There are only 4 floating point numbers on ^1 and 6 on ^2.


	We show ^d-1∩ P^d ⊈{-1,0,1}^d implies d ≥ 4.
	Suppose there is a non-trivial p ∈^d-1∩ P^d with d minimal.
	Let x_i/2^e_i denote the canceled fraction of its i-th coordinate.
	We have that all x_i ≠ 0, x_i are odd numbers and all e_i >0 (since p is not one of the 2d poles and d is minimal).
	
	W.l.o.g. e_1 ≤ e_2 ≤…≤ e_d.
	We rewrite the sphere equation 1 = ∑_j=1^d (x_i/2^e_i)^2 to
	
    x_1^2    = 4^e_1 - ∑_j=2^d4^e_1 - e_j x_j^2   .

	For an odd integer y, we have y^2=(2k+1)^2= 4(k^2+k)+1, leading to the congruence
	
    1 ≡ 0 - ∑_j=2^dχ_e_1(e_j)   4  .

	Where the characteristic function χ_e_1(e_j) is 1 for e_1 = e_j and 0 otherwise.
	For d∈{2,3} the right hand side can only have values of 0,-1 or -2, a contradiction.

Note that theorem <ref> translates to spheres with other radii through scaling.
Suppose a sphere in ^3 of radius 2^j has a non-trivial solution y ∈ P^3, then y/2^j ∈ P^3 and would be on ^2, too.







 §.§ Snapping to Rational Points

We now describe how to compute a good rational approximation exactly on the unit sphere ^d-1.
The input point x ∈^d can be given by any finite means that allows to compute rational approximations of arbitrary target precision – E.g. rational approximations of Cartesians for spherical coordinates.
For the input x, we denote its closest point on ^d-1 with x/‖ x ‖_2.
The stereographic projection τ and its inverse mapping σ provide σ(τ( x/‖ x‖_2 ) )= x/‖ x‖_2, since the argument is on ^d-1.
Instead of determining the value of τ exactly, we calculate an approximation y ∈^d and finally evaluate σ(y) under exact, rational arithmetics.
Hence, the result σ(y) is exactly on ^d-1.
r.35
	
	[scale=.13]
	(0,0) circle (10);
	
	(10,-10) circle (.25) node [right] 			x;
	(7,-7) circle (.25) node [left] 				x / ‖ x ‖_2;
	(8,-6) circle (.25) node [right] 				σ(y);
	(0,+10) circle (.25) node [above] 			(0,1);
	(4.14213562,0) circle (.25) node [above left] τ(x/‖ x ‖_2);
	(5,0) circle (.25) node [above right] 		y;
	
	[-] (-12,0) – node [above]  (12,0);	 
	[-,dashed,gray] (0,0) – node [above]  (10,-10); 
	[-] (0,10) – node [above]  (7.07106781,-7.07106781);
	[-] (0,10) – node [above]  (8,-6);
	

The stereographic projection does not preserve distances, leaving it open to bound the approximation error and the size of the resulting denominators.
We use the rotation symmetry of the sphere to limit the stretching of σ (c.f. Lemma <ref>):
For a non-zero point x ∈^d we can assume that i=d maximizes |x_i| and x_d < 0, otherwise we change the standard orthonormal basis by swapping dimension i and d and using a negative sign for dimension d.
Note that such rotations do not change the actual coordinate values.
To keep the size of denominators in σ(y) small, we use fixed-point arithmetics to determine y ∈^d-1 (c.f. Lemma <ref>).



See Algorithm <ref> for a precise description.
Note that the rational point y in statement <ref> solely needs to meet the target approximation in the individual coordinates for

    τ_i (x / ‖ x ‖_2 ) = x_i/‖ x ‖_2 - x_ d  .

Generally, this can be determined with methods of `approximate expression evaluation' to our target precision<cit.>.
If x is an approximation to a geo-referenced point, this denominator is well conditioned for calculations with multi-precision floating-point arithmetics<cit.>.
Using exact rational arithmetics for statement <ref>, we obtain a rational Cartesian coordinates on the unit sphere.

 
	For d > 1 and x ∈^d-1 with x_d = min_i -|x_i|, we have
	
    ‖τ(x) ‖_2 ≤√(√(d)-1/√(d)+1) < 1   .



Using x_d = min_i -|x_i| and ∑_ix_i^2= 1, we have the bounds 1/d ≤ x_d^2 ≤ 1 and

    ‖τ ( x) ‖_2^2
       = ∑_i=1^d-1 x_i^2/(1-x_d)^2
     = 1 - x_d^2/(1-x_d)^2
     = 1+x_d/1-x_d≤1-1/√(d)/1+1/√(d) .

Where the latter term is in (0,1) for any d.
Hence the (d-1)-ball ^d-1_1 = { x ∈^d  : ‖ x ‖_2 ≤ 1} contains τ(x).



 §.§ Approximation Quality

See <cit.> for an earlier version of this paper with a weaker, but elementary, analysis.


We consider the problem in  the 2D hyperplane H_pyy', defined by two points y=σ(x), y'=σ(x') on ^d-1 and the projection pole p∈^d.
Given the rotation step in Algorithm <ref>, the projection plane H_0={ x ∈^d: x_d = 0} separates p and y,y' in ^d and in H_pyy'.
Since each q ∈ H_0 ∩ S^d-1 has ‖ q - p‖_2=√(2) (consider pq in H_0pq), 
the circumcircle C of p,y and y' contains exactly two of these points.
Hence, the line of H_pyy'∩ H_0 is orthogonal to the circumcircle's diameter through p.
Moreover, the circles diameter is in [√(2),2].
We denote with x the point that is closer to p in H_pyy', meaning ‖ x‖_2 ≤‖ x'‖_2.
Note that x' and x can be on the same or opposite circumcircle halves.




	[scale=.3]
	(P) at (90:10);
	(V) at (210:10);
	(U) at (230:10);
	(M) at (0,0);
	(HP) at  (190:10);
	(HP2) at (-10:10);
	[gray,dashed] (0,0) circle (10);
	
	
	(P) circle (.25) node [above] p;
	(V) circle (.25) node [below left] y';
	(U) circle (.25) node [below left]  y;
	(M) circle (.25) node [right] m;
	
	[name path=PV] (P) – node [above]  (V);
	[name path=PU] (P) – node [above]  (U);
	[name path=XY] (V) – node [above]  (U);
	[name path=VM] (V) – (0,0) – (P);
	[name path=HPP] (HP) – (HP2)  node [below]  H_0 ;
	
	[name intersections=of=HPP and PU,by=HPP-PY];
	(V) – (U);
	[name path=XYSHIFT] (HPP-PY) – +((V)-(U));
	[name intersections=of=PV and XYSHIFT,by=PX-XYSHIFT];
	(PX-XYSHIFT) – (HPP-PY);
	[name intersections=of=PV and HPP,by=HPP-PX];
	(HPP-PX) – (HPP-PY);
	
	(HPP-PX) circle (.25) 				node [ below]x';
	(HPP-PY) circle (.25) 				node [below, xshift=6]x;
	((P)!(HPP-PY)!(V)) circle (.25) 	node [above left]		b;
	(PX-XYSHIFT) circle (.25) 			node [above, xshift=2]		a;
	
	[name path=B] ((P)!(HPP-PY)!(V)) – (HPP-PY);
	
	
	
	[dotted, gray] (P) – (HP)  node [midway,left ] =√(2);
	[dotted, gray] (P) – (HP2) node [midway,right] =√(2);
	
	[name path=MH0, dotted, gray] ((HP)!(M)!(HP2)) – (M);
	
	
	(0,2) arc (90:210:2);
	[above left] (M) γ;
	
	
	



In this section we denote with B=b x the perpendicular from x on py',E=xx', L=yy' and 
L_x=xa its triangle scaled version meeting x.
Note that B and L_x are above H_0, hence above E.

 
	For x,x' ∈^d-1_1 with ‖ x‖_2 ≤‖ x'‖_2, we have 
	
    ‖ x ‖_2/‖ p-σ(x) ‖_2‖σ(x) - σ(x') ‖_2 ≤‖ x - x' ‖_2 .




	We show L_x ≤ E by proofing α≤β for the two angles
	
    β   := ∡ b x x'
    α   := ∡ a x b .


	The inner angle sum of xab with a supplementary angle argument and triangle scaling provide ∡ py'y = 90^∘ + α.
	Let m denote the center of C.
	Since pm is orthogonal on H_0 and ∡ bx'x = 90^∘ - β, we have ∡ mp x' = β.
	In the isosceles triangle py'm, the central angle γ = 180^∘ - 2β.
	Fixing arc py' on C for the inscribed angle theorem provides ∡ y'yp=γ/2.
	
	Now, suppose α > β.
	The inner angle sum of pyy' states
	
    0 ≤∡ y'py
    	   = 180^∘ - ∡ y'yp - ∡ py'y
       = 180^∘ - ∡ y'yp - (90^∘ + α)
       = 180^∘ - γ/2  - (90^∘ + α)
        = -α + β

	a contradiction.



 
	For x, x' ∈^d-1_1 , we have
	
    ‖σ( x) - σ(x') ‖_2 ≤ 2  ‖ x - x' ‖_2.



	Using Lemma <ref>, we have L_x≤ E and the statement follows via triangle scaling:
	
    L = L_x  py / px≤ 2 L_x≤  2 E  ,

	since px≥ 1 and py≤ 2.

This statement is tight,
considering the two points x = 0 and
x'=( ε/√(d-1), …, ε/√(d-1)).
We have ‖ x - x'‖_2=ε and 
‖σ(x) - σ(x')‖_2
= 21/√(1+ε^2 )ε.


Algorithm <ref> calculates an ε-approximation exactly on the unit sphere.


Let x^* = x /‖ x ‖_2 and σ(y) denote the result.
Given the rotation, x^* holds for Observation <ref>.
Hence, we can use Lemma <ref> to derive

    ‖σ(y) - x^*‖_∞   = ‖σ(y) - σ(τ(x^*))‖_∞
       ≤‖σ(y) - σ(τ(x^*))‖_2 
       ≤ 2 ‖ y - τ(x^*) ‖_2 
       ≤2√((d-1) ε^2/4(d-1)) = ε

as upper bound on the approximation error.

This analysis is rather tight, as demonstrated by the red curve and points in Figure <ref>.




 §.§ Denominator Sizes

We now describe a relation between rational images of σ and the lowest common multiple of  denominators of its rational pre-images.
This leads to several strategies for achieving small denominators in the results of Algorithm <ref>.
 
Let x ∈^d-1∩^d-1_1 with x_i = p_i/q_i and Q = lcm(q_1, …, q_d-1) be the lowest common multiple, then 

    σ_k (x)= n_k/m
 
with integers n_i, m ∈{-2Q^2, … , 2Q^2} for all 1≤ k ≤ d.


Let q'_i ∈{1,…,Q} such that q'_i · q_i = Q for all i.
Since the formula of σ is similar in all but the last dimension, we describe the following two cases.
For k=d, we have

    σ_k (x)
       =-1 + ∑_i=1^d-1p_i^2 / q_i^2/1 + ∑_i=1^d-1p_i^2 / q_i^2
    =-Q^2 + ∑_i=1^d-1q'_i^2p_i^2 /Q^2 + ∑_i=1^d-1q'_i^2p_i^2 =: n_k/m

Using the bound x ∈^d-1_1, we have 0 ≤∑_i=1^d-1q'_i^2p_i^2 ≤ Q^2 and we derive for n_k and m

    |n_k|    = |-Q^2+ ∑_i=1^d-1q'_i^2p_i^2 | ≤ Q^2
    
    	m    = Q^2 + ∑_i=1^d-1q'_i^2p_i^2 ≤ 2Q^2

	
For k < d, we have

    σ_k (x)
    	   =2 p_k/q_k/1 + ∑_i=1^d-1p_i^2 / q_i^2
       = Q^2 · 2p_k/q_k /Q^2 + ∑_i=1^d-1q'_i^2 p_i^2 
    
    	   = Qq'_k · 2 p_k /Q^2 + ∑_i=1^d-1q'_i^2 p_i^2    =: n_k/m

Using the bound x ∈^d-1_1, we have that each |p_i| ≤ q_i and this bounds |n_k| = Qq'_k · 2 |p_k| ≤ 2Q^2.
We already discussed the bound on m in the first case.


Note that we apply this lemma in practice with fixed-point binary numbers p_i/q_i ∈ P_s.
Meaning all q_i = 2^s = Q for some significant size s.

	Denominators in ε-approximations of Algorithm <ref> are at most 
	
    10(d-1)/ε^2 .



Using standard multi-precision floating point arithmetics allows to derive rational values y, with denominators that are Q=⌈2√(d-1)/ε⌉.
Using ε≤ 1/8 and Lemma <ref> bounds the size of the denominators in images σ with

    2Q^2    ≤ 2(1+ 2√(d-1)/ε)^2
       = 2/ε^2( ε^2 +ε4√(d-1)_≤ (d-1) + 4(d-1) ) .



For certain dimensions and in practice(c.f. Section <ref>), we can improve on the simple usage of fixed-point binary numbers.
For ^1 we can rely on the continued fraction algorithm to derive rational approximations of α = τ(x/‖ x‖_2) with | α - p/q| ≤ 1/2q^2.
Using this in Algorithm <ref> leads to approximations with ε=1/q^2 on the circle ^1 with denominators of at most 2q^2.

Note that for ^d with d≥ 2 one can rely on algorithms for simultaneous Diophantine approximations (c.f. Theorem <ref>) to keep the lowest common multiple Q in Lemma <ref> small.
Note that it might well be simpler to find Diophantine approximations with small Q.

There have been many approaches to find generalizations of the continued fraction algorithm for d>1.
One of the first approaches is the Jacobi-Perron algorithm, which is rather simple to implement<cit.>(c.f. Section <ref>).
More advanced approaches <cit.> rely on the LLL-algorithm for lattice basis reduction<cit.>.
For d=2 there is an algorithm to compute all Dirichlet Approximations<cit.>, which we find hard to oversee given its extensive presentation.
Moreover, their experimental comparison shows that the Jacobi-Perron algorithm is practically well suited for d=2.

We close this section with a transfer result of Theorem <ref> with our Theorem <ref> and Lemma <ref>.

Let x ∈^d-1 and N ∈.
There is p ∈^d-1 and q ∈{1,…,N} with  

    ‖ x - σ(1/qp) ‖_∞≤2√(d-1)/q√(N)

and all denominators of σ(1/qp) are at most 2q^2.

This existence statement allows for brute-force computations.
However, we just use it for comparisons in Section <ref>.



§ IMPLEMENTATION

Apart from <cit.> for ^2, most implementations of spherical Delaunay triangulations are not `stable'.
Approaches based on d-dimensional convex hull algorithms produce only a tessellation for input not exactly on ^d-1. (c.f. Section <ref>)

Few available implementations allow dynamic point or constraint insertion and deletion – not even in the planar case of ^2.
The `Computational Geometry Algorithms Library' (CGAL <cit.>) is, to our knowledge, the sole implementation providing dynamic insertions/deletions of points and constraint line segments in ^2. 

With <cit.>, we provide open-source implementations of Algorithm <ref> for ^d. In <cit.>, we provide an implementation for spherical Delaunay triangulations on ^2 with ε-stable constructions of intersection points of constraint line-segments (c.f. Section <ref>).



 §.§ RATional Sphere Snapping for ^d

http://www.github.com/fmi-alg/libratssLibratss is a C++ library which implements Algorithm <ref>, based on the open-source GMP library for exact rational arithmetics <cit.> and the GNU `Multiple Precision Floating-Point Reliably'(MPFR) library<cit.>.
The implementation allows both, input of Cartesian coordinates of arbitrary dimension and spherical coordinates of ^2.
Note that this implementation allows geometric algorithms, as for d-dimen­sional convex hull, to rely on rational input points that are exactly on ^d-1.
In light of the discussion on the denominator sizes in Section <ref>, we provide two additional strategies to fixed-point snapping, as analyzed in Theorem <ref>.
We implemented the Continued Fraction Algorithm to derive rational ε-approximations with small denominators and the Jacobi-Perron algorithm for ^2.
The library interface also allows to automatically chose the approximation method which results in smaller denominators, approximation errors or other objectives, like byte-size.




 §.§ Incremental Constrained Delaunay Triangulation on ^2

http://www.github.com/fmi-alg/libdts2Libdts2 implements an adapter for the dynamic constraint Delaunay triangulation in the Euclidean plane ^2 of CGAL.
Since this implementation requires an initial outer face, we introduce an small triangle, that only contains the north-pole, to allow subsequent insertions of points and constraints.
For points exactly on the unit sphere, the predicate `is A in the circumcircle of B,C and D' reduces to the well studied predicate `is A above the plane through B,C and D'.
The implementation overloads all predicate functions accordingly and uses Algorithm <ref> for the construction of rational points on the sphere for intersections of Great Circle segments.



  §.§.§ ϵ-stable geometric constructions
 
Any means of geometric construction that allows to approximate a certain point, can be used as input for Algorithm <ref> – E.g. the intersection of Great Circle segments.
Consider two intersecting segments of rational points on ^2. 
The two planes, containing the segments and the origin as a third point, intersect in a straight line.
Each (rational) point on this line can be used as input for our method, as they identify the two intersection points on the sphere.
Using such input for Algorithm <ref> allows simple schemes to derive stable geometric constructions of rational points on ^d within a distance of ε to the target point.




§ EXPERIMENTS

We used real world and synthetic data for our experiments.
Geo-referenced data was sampled from regional extracts from the OpenStreetMap project<cit.>, as of January 26th, 2017.
Random Cartesian coordinates of points on ^d were created with the uniform generator 2 of <cit.>.
All benchmarks were conducted on a single core of an Intel Xeon E5-2650v4.
Peak memory usage and time were measured using the  utility.





 §.§ Approximation Quality and Size
 
We experimentally analyze the actual approximation error in results of Algorithm <ref> for several levels of ε using the MPFR library.
In this section e denotes the significands required in statement <ref> of Algorithm <ref> for the required result precision ε. This is 

    e = ⌈ -log_2 ( ϵ/2√(d-1)) ⌉  .

We simply setup the MPFR data types with significand sizes up to 1024 Bits, and conducted our experiments on much lower levels of e.
This allows us to derive some `measure' of the actual approximation errors of our method.

We analyzed the approximation errors δ and denominator bit-sizes q for 
100 random points on ^2. 
Figure <ref> compares the results of our algorithm under several levels of target precision e and strategies for statement <ref> in our method.
The magenta line indicates the quality and size of the approach in <cit.>.
The red line indicates the bounds of our Theorems <ref> and <ref> on the fixed-point strategy, while the yellow line indicates the bound of Corollary <ref>.
Note that results using the Jacobi-Perron strategy (blue dots) allows our method to further improve on the fixed-point strategy (red dots).
Note that we use Liouville's lower bound as statement on the approximability of a worst-case point.
There might well be points of higher algebraic degree that allow better approximations (c.f. Section <ref>).

Table <ref> exhibits average approximation errors δ, denominator bit-sizes q and the computation time t of our method for millions of points.
Synthetic data sets have several dimensions, while the real world data sets have dimension 3.
For ^2, we provide comparison of the fixed-point strategy (fx) with the Jacobi-Perron strategy (jp) of our method.
Using e=31 is sufficient to obtain results exactly on ^2 with a δ of less than 1cm, relative to a sphere with radius of the earth.
This is enough for most applications dealing with spatial data and allows storage within the word size of contemporary computing hardware.
This allows practical applications on ^2 to store 4 integer long values for the 3 numerators and the common denominator (c.f. Lemma <ref>) occupying 32 Bytes.
Note that storing 3 double values occupies 24 Bytes but cannot represent Cartesian coordinates exactly on the sphere.








 §.§ Constrained Delaunay Triangulation with Intersection Constructions
 
A Constrained Delaunay Triangulation of a point set contains required line-segments as edges, but is as close to the Delaunay triangulation as possible <cit.>.
We used very large street networks of several regions from the OpenStreetMap project for points and constraint edges – E.g. each line-segment of a street is an edge in the result triangulation.
Since ∼0.5% of the line-segments in these data sets intersect, we approximated the intersection points using e=31 for Algorithm <ref>.
Table <ref> exhibits total running time, peak memory usage and the result sizes of our  implementation.
Small data sets like Saarland and Germany allow quick calculation on a recent workstation computer.
See Figure <ref> for the Ecuador dataset.
Note that the current implementation has a storage overhead for each point, as we keep the results of the GMP library rather than truncating to integers of architectures word size.
Computing the triangulation for the planet data set was only possible on rather powerful hardware with at least 550 Gigabytes of memory taking half a day.



§ OPEN PROBLEMS

From a practical point of view, it is of great interest to bound the storage size of denominators to a maximum of 64Bits – the word size of current computing architectures.
We seek to improve our (already satisfactory) results by using advanced algorithms for simultaneous approximation, like the LLL-algorithm or the Dirichlet approximation algorithm for ^2.

For the theoretical part, we are interested if finding simultaneous rational approximations with small lowest common multiple of the denominators is simpler than finding Dirichlet approximations.
We are also interested in generalizing the method to provide rational approximations with small absolute errors on ellipsoids with rational semi-principal axes – E.g. the geographic WGS84 ellipsoid.

abbrv 
	
	
	

§ PROOF OF LIOUVILLE'S APPROXIMATION THEOREM
 
	This nice proof was translated from the German  https://de.wikibooks.org/wiki/Beweisarchiv:_Algebra:_K
		See <cit.> for the original proof in French language.
		
		equationsection
		
		
			Let α∈ be algebraic of degree n and root of the corresponding polynomial f(X)∈[X] of degree n, meaning
			
    f(α) = a_0 + a_1α + ⋯ +a_nα^n = 0
 with a_0, …, a_n ∈ and a_n 0.
			
			Polynomial division with the linear factor X-α in the ring [X] provides
			
    f(X) = (X-α)· g(X).

			
			Note that the polynomial g(X) has algebraic coefficients and is not necessarily in [X].
			However, the mapping →, t↦ g(t) is continuous, by means of real numbers c_1>0, c_2>0 with
			
    |g(x)| ≤ c_1

			for |α-x| < c_2.
			Since n<∞, we can assume w.l.o.g. that no additional roots are in this neighborhood of α, meaning
			
    f(x) 0

			for |α-x| < c_2 and xα.
			
			Claim: The statement of the Theorem holds for
			
    c:=min{c_2, 1/c_1}  .

			
			Suppose there are p,q∈, q>0 with
			
    |α-p/q | < c/q^n  .

			We show that his implies α= pq.
			
			From (<ref>), we immediately derive
			
    |α-p/q| < c ≤ c_2  ,

			leading (<ref>) to imply | g( pq)  | ≤ c_1.
			We derive, from (<ref>) and again (<ref>), that
			
    |f(p/q) | = |p/q -α|·|g(p/q)| < c/q^n· c_1 ≤1/q^n,

			meaning |q^n · f(p/q)| < 1.
			
			However q^n· f(p/q) = a_0q^n + a_1pq^n-1 + ⋯ + a_n p^n∈ and its absolute value is smaller than 1, hence has to be 0.
			Moreover, f( pq) = 0 and (<ref>) with (<ref>) imply α= pq, which closes the argument.
		
		
		
		

§ PROOF OF DIRICHLET'S APPROXIMATION THEOREM

		The folklore proof bases on Dirichlet's famous Pigeonhole Principle.
		See https://proofwiki.org/wiki/Dirichlet
			
				We consider the partition of [0,1]^d in N^d regular d-cubes of length L=√(N).
				We further define a sequence of points (a^(j))_j=1,…, N^d+1∈ [0,1]^d with a^(j) := j ·α - ⌊ j ·α⌋ (component wise operations).
				There are indices k > l such that the points a^(k) and a^(l) are contained in the same d-cube.
				We have the (component wise) inequalities
				
    -1/L < a^(k)   - a^(l) < 1/L
    
    				-1/L < kα - ⌊ k α⌋   - lα+⌊ l α⌋ < 1/L
    
    				-1/L < (k-l)α   - (⌊ k α⌋  - ⌊ l α⌋ ) < 1/L

				Setting q=k-l and p_i = ⌊ kα_i⌋ - ⌊ lα_i⌋ provides integers as required.
			
			

			

§ REDUCTIONS OF SPHERICAL PREDICATES TO CARTESIAN ORIENTATION PREDICATES

			We first describe a reduction from the spherical predicates to well studied Cartesian predicates.
			
			
				Let p_1, p_2∈𝕊^2 with p_1 ≠ p_2 and P the plane containing p_1,p_2 and the origin (0,0,0) and C be the Great Circle through p_1 and p_2. For q ∈𝕊^2 we have 
				
    q  left-of  P     q  left-of  C
    
    				q ∈ P     q ∈ C
    
    				q  right-of  P     q  right-of  C

			
			
				𝕊^2 ∩ P = C and 𝕊^2= L ∪ C ∪ R.
			
			
			
				Let P denote the plane through non-identical points p_1, p_2, p_3 ∈𝕊^2 and the half space containing the origin (0,0,0) is called `below P'.
				We further call S_123⊆^3 the closed volume of the sphere with p_1,p_2,p_3 and the origin on its surface. For a point q ∈𝕊^2 we have
				
    q  above  P      q ∈ S_123∖∂ S_123
    
    				q ∈ P 		        q ∈∂ S_123
    
    				q  below  P      q ∉ S_123

			
			
			
				P is uniquely determined because three different points on the unit sphere are not co-linear.
				Since S_123 and 𝕊^2 are spheres, their cuts with P are circles in P and the two circles are identical as they contain p_1,p_2 and p_3 on their boundary.
				This circle C has a radius of at most 1 and partitions the points of the unit sphere into three sets 
    𝕊^2 = A ∪ C ∪ B
 where A ⊆ S_123⊋ B.
				If C is a Great Circle we resolve ambiguity for `above' and the center of S_123 by choosing the open half spaces that first contain (0,0,1), then (0,1,0) and eventually (1,0,0).
				We have 𝕊^2 ≠∂ S_123, since the origin is a fourth point on S_123 and q ∈ P iff. q ∈ C iff. q ∈∂ S_123.
				Therefore it is sufficient to show q ∈ A	 q  above  P.
				To this end we consider the convex volume of the unit sphere S ⊆^3 and D = S ∩ S_123. Note that ∂ D contains C and A.
				Since the cut with the closed half-space of the plane P cuts a convex body into at most three parts and C ⊆ P, we have that all of A is `above' P.
			
			
			

		
