Algebraic Real Numbers ( real )

Definition

An instance x of the data type real is an algebraic real number. There are many ways to construct a real: either by conversion from double, bigfloat, integer or rational or by applying one of the arithmetic operators + , - ,*,/ or $\sqrt[d]{}$ to real numbers. One may test the sign of a real number or compare two real numbers by any of the comparison relations = ,! = , < , < = , > and > =. The outcome of such a test is exact. There is also a non-standard version of the sign function: the call x.sign(integer q) computes the sign of x under the precondition that | x| < = 2^-q implies x = 0. This version of the sign function allows the user to assist the data type in the computation of the sign of x, see the example below.

There are several functions to compute approximations of reals. The calls to_bigfloat(x) and x.get_bigfloat_error() return bigfloats xnum and xerr such that | xnum - x| < = xerr. The user may set a bound on xerr. More precisely, after the call x.improve_approximation_to(integer q) the data type guarantees xerr < = 2^-q. One can also ask for double approximations of a real number x. The calls to_double(x) and x.get_double_error() return doubles xnum and xerr such that | xnum - x| $\leq$ | xnum|*xerr. Note that xerr = $\infty$ is possible.

#include < LEDA/real.h >

Creation

reals may be constructed from data types double, bigfloat, long, int and integer. The default constructor real() initializes the real to zero.

Operations

double	to_double(real)	returns the current double approximation of x.
bigfloat	to_bigfloat(real)	returns the current bigfloat approximation of x.
double	x.get_double_error()	returns the relative error of the current double approximation of x, i.e., | x - to$\_$double(x)| is bounded by x.get$\_$double$\_$error()*| to$\_$double(x)|.
bigfloat	x.get_bigfloat_error()	returns the absolute error of the current bigfloat approximation of x, i.e., | x - to$\_$bigfloat(x)| is bounded by x.get_bigfloat_error().
int	sign(real x)	returns the sign of (the exact value of) x.
int	x.sign(integer q)	as above. Precondition if | x| < = 2-q then x = 0.
void	x.improve_approximation_to(integer q)
		(re-)computes the approximation of x such that x.get$\_$bigfloat$\_$error() < = 2-q after the call x.improve_approximation_to(q).
void	x.compute_with_precision(long k)
		(re-)computes the approximation of x; each numerical operation is carried out with k binary places.
void	x.guarantee_relative_error(long k)
		(re-)computes an approximation of x such that the relative error of the approximation is less than 2-k.
ostream&	ostream& O << x	writes the closest interval that is known to contain x to the output stream O
istream&	istream& I >> real& x	reads real number x from the output stream I (in double format)
real	sqrt(real x)	$\sqrt{x}$
real	root(real x, int d)	$\sqrt[d]{x}$ Precondition d > = 2
real	abs(real x)	absolute value of x
real	sqr(real x)	square of x
real	dist(real x, real y)	euclidean distance of vector (x,y) to the origin
real	powi(real x, int n)	n.th power of x

Implementation

A real is represented by the expression which defines it and a double approximation x' together with a relative error bound $\epsilon_{x}^{}$. The arithmetic operators +, -, *, /, $\sqrt[d]{}$ take constant time. When the sign of a real number needs to be determined, the data type first computes a number q, if not already given as an argument to sign, such that | x| < = 2^-q implies x = 0. The bound q is computed as described in [55].

Using bigfloat artihmetic, the data type then computes an interval I of maximal length 2^-q that contains x. If I contains zero, then x itself is equal to zero. Otherwise, the sign of any point in I is returned as the sign of x.

Two shortcuts are used to speed up the computation of the sign. Firstly, if the double approximation already suffices to determine the sign, then no bigfloat approximation is computed at all. Secondly, the bigfloat approximation is first computed only with small precision. The precision is then doubled until either the sign can be decided (i.e., if the current approximation interval does not contain zero) or the full precision 2^-q is reached.

Example

We give two examples of the use of the data type real. The examples deal with the Voronoi diagram of line segments and the intersection of line segments, respectively.

The following incircle test arises in the computation of Voronoi diagrams of line segments [15,13]. For i, 1 < = i < = 3, let l_i : a_ix + b_iy + c_i = 0 be a line in two-dimensional space and let p = (0, 0) be the origin. In general, there are two circles passing through p and touching l₁ and l₂. The centers of these circles have homogeneuos coordinates (x_v, y_v, z_v), where

$$\begin{eqnarray*} x_v & = & a_1 c_2 +a_2 c_1 \pm\mbox{sign}(s)\sqrt{2 c_1 c_2 (\... ... (\sqrt{N} - D)} \\ [2mm] z_v & = & \sqrt{N} - a_1 a_2 - b_1 b_2 \end{eqnarray*}$$

and

$$\begin{array}{cccccc} s & = & b_1 D_2 - b_2 D_1, & N & = & (... ...= & a_1 D_2 - a_2 D_1, & D & = & a_1 a_2 - b_1 b_2. \end{array}$$

Let us concentrate on one of these (say, we take the plus sign in both cases). The test whether l₃ intersects, touches or misses the circle amounts to determining the sign of

$$\begin{eqnarray*} E := dist^{2}(v,l_3) - dist^{2}(v,p) = \frac{(a_3 x_v + b_3 y_v + c_3)^2}{a_3^2 + b_3^2} - (x_v^2 + y_v^2). \end{eqnarray*}$$

The following program computes the sign of $\tilde{E}$ : = (a₃² + b₃²)*E using our data type real.

int INCIRCLE( real a₁, real b₁, real c₁, real a₂, real b₂, real c₂, real a₃, real b₃, real c₃ )
{
real RN = sqrt((a₁*a₁ + b₁*b₁)*(a₂*a₂ + b₂*b₂));
real RN₁ = sqrt(a₁*a₁ + b₁*b₁);
real RN₂ = sqrt(a₂*a₂ + b₂*b₂);
real A = a₁*c₂ + a₂*c₁;
real B = b₁*c₂ + b₂*c₁;
real C = 2*c₁*c₂;
real D = a₁*a₂ - b₁*b₂;
real s = b₁*RN₂ - b₂*RN₁;
real r = a₁*RN₂ - a₂*RN₁;
int sign_x = sign(s);
int sign_y = sign(r);
real x_v = A + sign_x* sqrt(C*(RN + D));
real y_v = B - sign_y* sqrt(C*(RN - D));
real z_v = RN - (a₁*a₂ + b₁*b₂);
real P = a₃*x_v + b₃*y_v + c₃*z_v;
real D₃² = a₃*a₃ + b₃*b₃;
real R² = x_v*x_v + y_v*y_v;
real E = P*P - D₃²*R²;
return sign(E);
}

We can make the above program more efficient if all coefficients a_i, b_i and c_i, 1 < = i < = 3, are k bit integers, i.e., integers whose absolute value is bounded by 2^k - 1. In [15,13] we showed that for $\tilde{E}$! = 0 we have |$\tilde{E}$| > = 2^{-24k - 26}. Hence we may add a parameter int k in the above program and replace the last line by

${\bf return} \; E.\mbox{sign}(24*k+26).$
Without this assistance, reals automatically compute a weaker bound of |$\tilde{E}$| > = 2^{-56k - 161} for $\tilde{E}$! = 0 by [14].

We turn to the line segment intersection problem next. Assume that all endpoints have k-bit integer homogeneous coordinates. This implies that the intersection points have homogeneous coordinates (X, Y, W) where X, Y and W are (4 k + 3) - bit integers. The Bentley-Ottmann plane sweep algorithm for segment intersection [59] needs to sort points by their x-coordinates, i.e., to compare fractions X₁/W₁ and X₂/W₂ where X₁, X₂, W₁, W₂ are as above. This is tantamount to determining the sign of the 8k + 7 bit integer X₁*W₂ - X₂*W₁. If all variables X_i, W_i are declared real then their sign test will be performed quite efficiently. First, a double approximation is computed and then, if necessary, bigfloat approximations of increasing precision. In many cases, the double approximation already determines the sign. In this way, the user of the data type real gets the efficiency of a floating point filter [31,60] without any work on his side. This is in marked contrast to [31,60] and will be incorporated into [59].