The data type real ( 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 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| < = 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()	returns the current double approximation of x.
double	to_float()	as above.
bigfloat	to_bigfloat()	returns the current bigfloat approximation of x.
double	x.get_double_error()	returns the absolute error of the current double approximation of x, i.e., | x - to$\_$double(x)| < = x.get$\_$double$\_$error().
bigfloat	x.get_bigfloat_error()	returns the absolute error of the current bigfloat approximation of x, i.e., | x - to$\_$bigfloat(x)| < = 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: The user guarantees that | x| < = 2-q is only possible if x = 0. This advanced version of the sign function should be applied only by the experienced user. It gives an improvement over the plain sign function only in some cases.
void	x.improve_approximation_to(integer q)
		recomputes the approximation of x if necessary; the resulting error of the bigfloat approximation satisfies x.get$\_$bigfloat$\_$error() < = 2-q
void	x.compute_with_precision(long k)
		recomputes the bigfloat approximation of x, if necessary; each numerical operation is carried out with a mantissa length of k. Note that here the size of the resulting x.get_bigfloat_error() cannot be predicted in general.
void	x.guarantee_relative_error(long k)
		recomputes an approximation of x, if necessary; the relative error of the resulting bigfloat approximation is less than 2-k, i.e., x.get$\_$bigfloat$\_$error() < = | x|*2-k.
ostream&	ostream& O << x	writes a bigfloat approximation of x to the output stream O. Note that the exact representation of the x instance is lost in the stream output and only the precision set by real::io_dec_precision is guaranteed.
istream&	istream& I >> real& x	reads x number x from the output stream I in double format. Note that stream input is currently impossible for more general types of reals.
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 an interval inclusion I that contains the exact value of the realThe 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 [,]. 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 initial interval 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 roughly 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. This procedure makes the sign computation of a real number x adaptive in the sense that the running time of the sign computation depends on the complexity of x.

Example

We give two typical examples for the use of the data type real that arise in Computational geometry. We admit that a certain knowledge about Computational geometry is required for their full understanding. The examples deal with the Voronoi diagram of line segments and the intersection of line segments, respectively.

The following incircle test is used 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 (... ...2 (\sqrt{N} - D)} \\ 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^{2} + b_1^{... ...2}) \\ r & = & 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 = \mbox{\small sqrt}((a_1*a_1+b_1*b_1)*(a_2*a_2+b_2*b_2))$;
real $RN_1 = \mbox{\small sqrt}(a_1*a_1+b_1*b_1);$
real $RN_2 = \mbox{\small sqrt}(a_2*a_2+b_2*b_2);$
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 * \mbox{\small sqrt}(C*(RN+D));$
real $y_v = B - sign_y * \mbox{\small 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 [].

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 boils down 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, an interval approximation is computed and then, if necessary, bigfloat approximations of increasing precision. In many cases, the interval approximation already determines the sign. In this way, the user of the data type real gets nearly the efficiency of a hand-coded 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].