6.4. Algebraic real numbers (class real)

[Warning]Warning

This chapter, section, or subsection is “in statu nascendi”. Therefore, it may be incomplete, or may appear in a different place in the final version, or may vanish completely.

Sometimes, even the arbitrarily high precision of bigfloats and the infinitely high precision of rationals is not enough. TODO say something about exact arithmetic, give motivation.

As we know from elementary number theory, the numbers on the real axis can be divided into the integer numbers, the rational numbers, the algebraic numbers, and the transcedental numbers. Let us briefly review this classification: The integer numbers are the numbers 0, 1, 2, 3, ..., and their negative counterparts -1, -2, -3, ... The rational numbers are the quotients p/q of two integer numbers p and q. We already know LEDA's classes integer for arbitrarily long integer numbers and rational and bigfloat for rational numbers. Each integer number is a rational number, but not each number on the real axis is rational:

An algebraic number is the root of a polynomial with rational coefficients. For example, the Golden Ratio Φ is a root of the polynomial Φ2 - Φ - 1. Each rational number is algebraic, but the converse is not true: Φ, for example, is not rational; it is thus a real algebraic number.

Algebraic numbers show up mostly in geometry. For example, the diagonal in a square whose sides are of length 1 has the length √2 (which is not rational either). To carry out exact computations with algebraic numbers, which mostly has to be done in computational geometry, LEDA provides the class real.

The representation of an algebraic number as the root of polynomials is “more complicated” than the representation of a rational number as the quotient of two integers. It is not the “most complicated” type of real numbers, though: The transcendental numbers are even more complicated than the algebraic numbers - they do not even have a representation as the root of a polynomial. For example, the circle constant π and Euler's number e are transcendental numbers. LEDA cannot deal with transcendental numbers, at least not without loss of precision - there is no number type class in LEDA that could represent π or e exactly.

A problem in which doubles are not enough

Let us consider the problem of writing a procedure that decides whether 3 circles intersect in a common point. Let C((a,b),r) be the circle with center point (a,b) and radius r, where all of a, b, and r are rational numbers. For example, the circles C((-4,0),5), C((4,0),5), and C((0,6),3) all intersect in the point (0,3), see (TODO there will be a drawing here making clear that (0,3) is an intersection point and depicting all variables of the algorithm below.)

Why are real algebraic numbers involved in this computation? The circle C((a,b),r) is the set of all points (x,y) that satisfy the equation

(x-a)2 + (y-b)2 = r2.

Therefore, the intersection of two circles is the solution of a system of two polynomial equations of degree 2, that is, the solution of a quadratic equation. This means that, in general, the solution will contain square roots that cannot be eliminated.

The following is a pseudo-code algorithm for computing the intersection of two circles C((a,b),r) and C((c,d),s)[45]:

Input a, b, r, c, d, s

e = c - a
f = d - b
g = sqrt(e2 + f2)
h = (g2 + r2 - s2)/(2g)

x1 = a + eh/g - (f/g)·sqrt(r2 - h2)
y1 = b + fh/g + (e/g)·sqrt(r2 - h2)

x2 = a + eh/g + (f/g)·sqrt(r2 - h2)
y2 = b + fh/g - (e/g)·sqrt(r2 - h2)

Output (x1,y1), (x2,y2)

Again, see TODO for an illustration of the quantities showing up in this algorithm.

Now, to decide whether three circles intersect in a common point, we intersect the circles pairwise and look for a point that occurs in each intersection.

In a trial implementation with C++'s floating point type double, we assume the “normal” case that we get two intersection points for each intersection, that is, we ignore the cases that there might be no intersection point, exactly one intersection point, and congruence:

Filename: CircleCircleIntersectionWithDoubles.C
LEDA users as of version 5.0: Note that header files are now partitioned into modules (subfolders)!
#include <cstdio>
#include <cmath>

// Computes the intersections of C((a,b),r) with C((c,d),s)

void CircleCircleIntersection( double a, double b, double r, 
                               double c, double d, double s )
{
  double e, f, g, h;
  double x1, y1, x2, y2;

  e = c - a; 
  f = d - b;                          
  g = sqrt(e*e + f*f); 
  h = (g*g + r*r - s*s) / (2*g); 
                              
  x1 = a + e*h/g + (f/g) * sqrt(r*r - h*h);
  y1 = b + f*h/g - (e/g) * sqrt(r*r - h*h);

  x2 = a + e*h/g - (f/g) * sqrt(r*r - h*h);
  y2 = b + f*h/g + (e/g) * sqrt(r*r - h*h);
  
  printf("(x1,y1) = (%.20lf,%.20lf)\n", x1, y1);
  printf("(x2,y2) = (%.20lf,%.20lf)\n", x2, y2);
}


int main()
{
  CircleCircleIntersection(-4, 0, 5, 4, 0, 5);
  CircleCircleIntersection( 4, 0, 5, 0, 6, 3);
  CircleCircleIntersection(-4, 0, 5, 0, 6, 3);
}

The output of this program on the author's machine is:

(x1,y1) = (0.00000000000000000000,-3.00000000000000000000)
(x2,y2) = (0.00000000000000000000,3.00000000000000000000)
(x1,y1) = (2.76923076923076960654,4.84615384615384581224)
(x2,y2) = (-0.00000000000000044409,2.99999999999999955591)
(x1,y1) = (0.00000000000000044409,2.99999999999999955591)
(x2,y2) = (-2.76923076923076960654,4.84615384615384581224)

By algebraic reasoning or by geometric insight from TODO, we know that three of these points should be equal to (0,3). Obviously, the program computes six different points. The differences in the outcome are caused by rounding errors introduced in the partial steps of the computation.

Even though three of the computed points are very close to one another, this closeness is of no value here: The problem was to decide whether three of them are equal, and a program based on doubles will say “no”. This is what exact arithmetic is all about: “Nearly equal” is the same as “completely different”. If two quantities are recognized as “nearly equal” although they should be “completely equal”, the computation is worthless.

The conclusion that we have to draw from this trial implementation is that doubles are not suitable for this intersection test.

The class real

LEDA can help us here: Implementation with real (ignoring the cases of no intersection point, exactly 1 intersection point, and congruence):

Filename: ThreeCirclesIntersection.C
LEDA users as of version 5.0: Note that header files are now partitioned into modules (subfolders)!
#include <LEDA/real.h>
#include <LEDA/list.h>
#include <iostream>

using leda::real;
using leda::list;
using std::cout;

struct point { real x, y; };

// Computes the intersections of C((a,b),r) with C((c,d),s)

list<point> CircleCircleIntersection( real a, real b, real r, 
                                      real c, real d, real s )
{
  real e, f, g, h;
  real x1, y1, x2, y2;
  list<point> result;
  point pint; // Point of intersection, if any

  e = c - a; 
  f = d - b;           
  g = sqrt(e*e + f*f); 
  h = (g*g + r*r - s*s) / (2*g); // Exercise: What if g == 0 ?
  
  if( r*r - h*h < 0 ) return result; // empty list: no intersection

  x1 = a + e*h/g + (f/g) * sqrt(r*r - h*h);
  y1 = b + f*h/g - (e/g) * sqrt(r*r - h*h);

  pint.x = x1; 
  pint.y = y1;
  result.push( pint );

  // r*r - h*h == 0 iff the circles touch each other in 1 point
  // In this case, we can return now; result list has only 1 element

  if( r*r - h*h == 0 ) return result;
                            
  x2 = a + e*h/g - (f/g) * sqrt(r*r - h*h);
  y2 = b + f*h/g + (e/g) * sqrt(r*r - h*h);

  pint.x = x2; 
  pint.y = y2;
  result.push( pint );

  return result;
}


// Tests whether 3 circles C((a,b),r), C((c,d),s) and C((e,f),t)
// have a common intersection point. If so, returns this point in pint

bool ThreeCirclesIntersection( real a, real b, real r, 
			       real c, real d, real s,
			       real e, real f, real t,
			       point& pint )
{
  list<point> L = CircleCircleIntersection( a, b, r, c, d, s );

  if( L.empty() ) return false;

  list<point> M = CircleCircleIntersection( a, b, r, e, f, t );

  if( M.empty() ) return false;

  // Exercise: What to do if L or M or both contain 1 point only?

  point l1 = L.head();
  point l2 = L.back();
  point m1 = M.head();
  point m2 = M.back();

  if( (l1.x == m1.x && l1.y == m1.y) 
      || (l1.x == m2.x && l1.y == m2.y) ) {
    pint.x = l1.x;
    pint.y = l1.y;
    return true;
  }

  if( (l2.x == m1.x && l2.y == m1.y)
      || (l2.x == m2.x && l2.y == m2.y) ) {
    pint.x = l2.x;
    pint.y = l2.y;
    return true;
  }

  return false;
}


int main()
{
  point pint;

  if( ThreeCirclesIntersection( -4, 0, 5, 4, 0, 5, 0, 6, 3, pint ) ) {
    cout << "Intersection point at ";
    cout << "(" << pint.x.to_double() 
	 << "," << pint.y.to_double() << ")\n";
  }
  else {
    cout << "No intersection point.\n";
  }
}

Output is

Intersection point at (0,3)

TODO In obigem Program gefällt mir die Definition von struct point nicht (sowas gibt's schon in LEDA, warum nicht nehmen?), und für die Schnittmengenberechnung würde man eine leda::set nehmen (aber dann muss der Parameter-Typ linear geordnet sein, was das Programm aufbläht; der zusätzliche Code könnte das, worum es geht, verschleiern.)

(Here will be the description of real's interface.)

Say something about the main usage: Used in programs (mainly from geometry) that make decisions based on the sign of an expression. Not used mainly for computing a number x with an accuracy of 5,000 digits after the comma. But: If x is needed in this accuracy to make a comparison, then it is automatically (and correctly!) computed.

TODO What to say about outputting a real as a decimal value? I'm afraid that I have to say something about the implementation (because of the term "current approximation" that shows up in the manual page description of to_double()).

reals vs. rationals and doubles

An attentive reader will bring forward the argument that rationals are sufficient in the above program: If three circles with rational coordinates and rational radius intersect in one point, then this point is rational. (Why? It is the interscetion of the perpendicular bisectors!)

So it could be solved by solving a system of rational linear equations, using the fact that that every double is a rational. But: The formulae are complex and this is laborious to code. And see Arno's argumentation below.

The above algorithm could also be implemented using only rationals: Use repeated squaring to successively eliminate all squares in the above formulae. But:

Diskussion mit Arno: Warum reals und nicht rationals? Die numerisch einfachen Fälle sind billig auszurechnen im Vgl. zu Repeated Squaring (exponentieller Aufwand in der Anzahl der Wurzeln des Ausdrucks) oder symbolischem Rechnen. "Numerisch einfach" bedeutet: Gerechnet wird bei reals eigentlich nur bei Vergleichen. (Alle anderen Operationen bestehen mehr oder weniger aus dem Aneinanderhängen von Bäumen.) Beide Bäume werden ausgewertet mit zu erzielender Endgenauigkeit. Wenn Vergleich dann schon mit Sicherheit ein korrektes Ergebnis liefert: ENDE. Ansonsten: Rechne nochmals mit höherer Endgenauigkeit usw., bis Vergleich Ergebnis liefert. "Numerisch einfach" = schon rel. geringe Endgenauigkeit genügt für Vergleich.

Das stimmt umso mehr, je mehr ein Algorithmus ausgehend von einer einmal berechneten Größe (die eine Wurzel enthält) weitere Größen berechnet und somit immer mehr Wurzeln ineinander schachtelt. Repeated Squaring ist dann nicht mehr praktikabel.

Warum oben (und auch sonst) nicht einfach "double < eps?" testen und nur mit doubles arbeiten? Problem dabei: Was ist das richtige epsilon? Es ist nicht a-priori für alle Eingaben zu entscheiden.

Summa summarum: reals sind exakt, bequem und eventuell sogar effizienter als rational. Doubles sind für exaktes Rechnen unbrauchbar.

Exercises

Exercise 82.

Extend the above program such that it also handles the cases that two or all three circles are congruent, or have the same center, or have exactly one intersection point.

Exercise 83.

TODO Represent Φ with a real.



[45] The proof is tedious and quite much math. Let

I)  (x-a)2 + (y-b)2 = r2
II) (x-c)2 + (y-d)2 = s2

be the circles. Just substituting I) into II) and then solving the resulting quadratic equation becomes really ugly. Here is a better way:

First, we translate the system such that the center (a,c) of the first circle is mapped to the origin (0,0): x → x + a, y → y + a. Thus we get the simplified system

I) x2 + y2= r2
II) (x-e)2 + (y-f)2 = s2

where e := c - a and f := d - b is the difference of the circles' x-coordinates and y-coordinates, respectively.

Next, we rotate by an angle -φ such that the center of the second circle, now positioned at (e,f), is mapped to (0,g) with g = sqrt(e2 + f 2). Note that g is the distance of the two centers and that cos(φ)=e/g and sin(φ)=f/g, which will come in handy later.

We are now left with the system

I) x2 + y2 = r2
II) (x-g)2 + y2 = s2

which we solve for x: By expanding II), substituting x2+y2 in I), and solving for x, we get x = -+ (g2+r2-s2)/(2g). Now substituting x in I) and solving for y, we see that y = +- 1/(2g) (r2 - (g2+r2-s2/(2g))2 ).

To simplify the formulae, we define h := g2+r2-s2/(2g). Note that h is the distance of center 1 and the line that connects the intersection points.

In a final step, we rotate back by φ and translate back by (a,b). The translation matrix is

| cos(φ)  -sin(φ)|
|                  |
| sin(φ)  cos(φ)|

with cos(φ)=e/g and sin(φ)=f/g. So we finally get the coordinates

x = a + he/g -+ sqrt(r2 - h2)
y = b + hf/g +- sqrt(r2 - h2)

of the intersection points.