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 `bigfloat`

s
and the infinitely high precision of `rational`

s 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.

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} = r^{2}.

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(e^{2}+ f^{2}) h = (g^{2}+ r^{2}- s^{2})/(2g) x1 = a + eh/g - (f/g)·sqrt(r^{2}- h^{2}) y1 = b + fh/g + (e/g)·sqrt(r^{2}- h^{2}) x2 = a + eh/g + (f/g)·sqrt(r^{2}- h^{2}) y2 = b + fh/g - (e/g)·sqrt(r^{2}- h^{2}) 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:

#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 `double`

s 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.

LEDA can help us here: Implementation with
`real`

(ignoring the cases of no intersection
point, exactly 1 intersection point, and congruence):

#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()`

).

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.

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

I) (x-a)^{2} + (y-b)^{2} = r^{2}

II) (x-c)^{2} + (y-d)^{2} = s^{2}

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) x^{2} + y^{2}= r^{2}

II) (x-e)^{2} + (y-f)^{2} = s^{2}

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(e^{2}
+ 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) x^{2} + y^{2} = r^{2}

II) (x-g)^{2} + y^{2} = s^{2}

which we solve for x: By expanding II), substituting
x^{2}+y^{2} in I), and
solving for x, we get x = -+
(g^{2}+r^{2}-s^{2})/(2g). Now substituting
x in I) and solving for y, we see that y = +- 1/(2g) (r^{2} - (g^{2}+r^{2}-s^{2}/(2g))^{2} ).

To simplify the formulae, we define h :=
g^{2}+r^{2}-s^{2}/(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(r^{2} - h^{2})

y = b + hf/g +- sqrt(r^{2} - h^{2})

of the intersection points.