Definition
An instance x of the data type real is a real algebraic number. There are many ways to construct a real: either by conversion from double, bigfloat, integer or rational, by applying one of the arithmetic operators + , - ,*,/ or to real numbers or by using the -operator to define a real root of a polynomial over 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 mathematically exact. We give consider an example expression to clarify this:
There are several functions to compute approximations of reals. The calls x.to_bigfloat() 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 x.to_double() and x.get_double_error() return doubles xnum and xerr such that | xnum - x| < = xerr. Note that xerr = is possible.
#include < LEDA/numbers/real.h >
Types
typedef polynomial<real> Polynomial | the polynomial type. |
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 | x.to_double() | returns the current double approximation of x. |
double | x.to_double(double& error) | |
as above, but also computes a bound on the absolute error. | ||
bigfloat | x.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 - x.todouble()| < = x.getdoubleerror(). |
bigfloat | x.get_bigfloat_error() | returns the absolute error of the current bigfloat approximation of x, i.e., | x - x.tobigfloat()| < = x.getbigfloaterror(). |
bigfloat | x.get_lower_bound() | returns the lower bound of the current interval approximation of x. |
bigfloat | x.get_upper_bound() | returns the upper bound of the current interval approximation of x. |
rational | x.high() | returns a rational upper bound of the current interval approximation of x. |
rational | x.low() | returns a rational lower bound of the current interval approximation of x. |
double | x.get_double_lower_bound() | |
returns a double lower bound of x. | ||
double | x.get_double_upper_bound() | |
returns a double upper bound of x. | ||
bool | x.possible_zero() | returns true if 0 is in the current interval approximation of x |
integer | x.separation_bound() | returns the separation bound of x. |
integer | x.sep_bfmss() | returns the k-ary BFMSS separation bound of x. |
integer | x.sep_degree_measure() | returns the degree measure separation bound of x. |
integer | x.sep_li_yap() | returns the Li / Yap separation bound of x. |
void | x.print_separation_bounds() | |
prints the different separation bounds of x. | ||
bool | x.is_general() | returns true if the expression defining x contains a -operator, false otherwise. |
bool | x.is_rational() | returns true if the expression is rational, false otherwise. |
rational | x.to_rational() | returns the rational number given by the expression. Precondition is_rational() has is true. |
int | x.compare(const real& y) | returns the sign of x-y. |
int | compare_all(const growing_array<real>& R, int& j) | |
compares all elements in R. It returns i such that R[i] = R[j] and i! = j. Precondition: Only two of the elements in R are equal. [ Experimental ] | ||
int | x.sign() | returns the sign of (the exact value of) x. |
int | x.sign(const 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(const integer& q) | |
recomputes the approximation of x if necessary; the resulting error of the bigfloat approximation satisfies x.getbigfloaterror() < = 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.getbigfloaterror() < = | x|*2^{-k}. | ||
ostream& | ostream& O « const real& x | |
writes the closest interval that is known to contain x to the output stream O. Note that the exact representation of x is lost in the stream output. | ||
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(const real& x) | |
real | root(const real& x, int d) | |
, precondition: d > = 2 |
Note: The functions real_roots and diamond below are all experimental if they are applied to a polynomial which is not square-free.
int | real_roots(const Polynomial& P, list<real>& roots, algorithm_type algorithm, bool is_squarefree) | |
returns all real roots of the polynomial P. | ||
int | real_roots(const Polynomial& P, growing_array<real>& roots, algorithm_type algorithm, bool is_squarefree) | |
same as above. | ||
int | real_roots(const int_Polynomial& iP, list<real>& roots, algorithm_type algorithm = isolating_algorithm, bool is_squarefree = true) | |
returns all real roots of the polynomial iP. | ||
real | diamond(int j, const Polynomial& P, algorithm_type algorithm, bool is_squarefree) | |
returns the j-th smallest real root of the polynomial P. | ||
real | diamond(rational l, rational u, const Polynomial& P, algorithm_type algorithm, bool is_squarefree) | |
returns the real root of the polynomial P which is in the isolating interval [l,u]. | ||
real | diamond_short(rational l, rational u, const Polynomial& P, algorithm_type algorithm, bool is_squarefree) | |
returns the real root of the polynomial P which is in the
isolating interval [l,u].
Precondition (u - l ) < 1/4 |
||
real | diamond(int j, const int_Polynomial& iP, algorithm_type algorithm = isolating_algorithm, bool is_squarefree = true) | |
returns the j-th smallest real root of the polynomial iP. | ||
real | diamond(rational l, rational u, const int_Polynomial& iP, algorithm_type algorithm = isolating_algorithm, bool is_squarefree = true) | |
returns the real root of the polynomial iP which is in the isolating interval [l,u]. | ||
real | abs(const real& x) | absolute value of x |
real | sqr(const real& x) | square of x |
real | dist(const real& x, const real& y) | |
euclidean distance of point (x,y) to the origin | ||
real | powi(const real& x, int n) | |
x^{n}, i.e., n.th power of x | ||
integer | floor(const real& x) | returns the largest integer smaller than or equal to x. |
integer | ceil(const real& x) | returns the smallest integer greater than or equal to x. |
rational | small_rational_between(const real& x, const real& y) | |
returns a rational number between x and y with the smallest available denominator. Note that the denominator does not need to be strictly minimal over all possible rationals. | ||
rational | small_rational_near(const real& x, double eps) | |
returns small_rational_between(x-eps, x+eps). |
Implementation
A real is represented by the expression which defines it and an interval inclusion I that contains the exact value of the real. The arithmetic operators + , - ,*, 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 [79]. Using bigfloat arithmetic, 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 [17,14]. For i, 1 < = i < = 3, let l_{i} : a_{i}x + b_{i}y + 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_{1} and l_{2}. The centers of these circles have homogeneuos coordinates (x_{v}, y_{v}, z_{v}), where
int INCIRCLE( real a_{1}, real b_{1}, real c_{1}, real a_{2}, real b_{2}, real c_{2}, real a_{3}, real b_{3}, real c_{3} )
{
real RN = sqrt((a_{1}*a_{1} + b_{1}*b_{1})*(a_{2}*a_{2} + b_{2}*b_{2}));
real RN_{1} = sqrt(a_{1}*a_{1} + b_{1}*b_{1});
real RN_{2} = sqrt(a_{2}*a_{2} + b_{2}*b_{2});
real A = a_{1}*c_{2} + a_{2}*c_{1};
real B = b_{1}*c_{2} + b_{2}*c_{1};
real C = 2*c_{1}*c_{2};
real D = a_{1}*a_{2} - b_{1}*b_{2};
real s = b_{1}*RN_{2} - b_{2}*RN_{1};
real r = a_{1}*RN_{2} - a_{2}*RN_{1};
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_{1}*a_{2} + b_{1}*b_{2});
real P = a_{3}*x_{v} + b_{3}*y_{v} + c_{3}*z_{v};
real D_{3}^{2} = a_{3}*a_{3} + b_{3}*b_{3};
real R^{2} = x_{v}*x_{v} + y_{v}*y_{v};
real E = P*P - D_{3}^{2}*R^{2};
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 [17,14] we showed that for ! = 0 we have || > = 2^{-24k-26}. Hence we may add a parameter int k in the above program and replace the last line 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 [65] needs to sort points by their x-coordinates, i.e., to compare fractions X_{1}/W_{1} and X_{2}/W_{2} where X_{1}, X_{2}, W_{1}, W_{2} are as above. This boils down to determining the sign of the 8k + 7 bit integer X_{1}*W_{2} - X_{2}*W_{1}. 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 [35,66] without any work on his side. This is in marked contrast to [35,66] and will be incorporated into [65].