The Dangers of Floating Point Arithmetic

The following two examples illustrate the dangers of floating point arithmetic in geometric computation.

Convex Hulls

This example was suggested by Stefan Schirra. We define a segment s and construct a set L of points consisting of the endpoints of s and the intersections between s and some random number of lines.

Since all points in L lie on s, the convex hull of L should have exactly two vertices. However, a sample run generates more than two vertices. The explanation is that, when the intersection between s and a line l is computed with the floating point kernel, the point of intersection does not necessarily lie on s but only near s.

#include <LEDA/geo/point.h>
#include <LEDA/geo/segment.h>
#include <LEDA/core/list.h>
#include <LEDA/geo/line.h>
#include <LEDA/geo/geo_alg.h>

using namespace leda;

int main()
{
  point p0(-LEDA_PI,-LEDA_PI);
  point p1(+LEDA_PI,+LEDA_PI);
  std::cout << "p0: (" << p0.xcoord() << "," << p0.ycoord() << ")\n";
  std::cout << "p1: (" << p1.xcoord() << "," << p1.ycoord() << ")\n";

  segment s(p0,p1);

  list<point> L; L.append(p0); L.append(p1);

  
  for (int i=0; i<10000;i++) {
    double ax,ay;
    rand_int >> ax; rand_int >> ay; point p(ax*LEDA_PI,ay*LEDA_PI);
    rand_int >> ax; rand_int >> ay; point q(ax*LEDA_PI,ay*LEDA_PI);
    line l(p,q); point r;
    if (l.intersection(s,r)) L.append(r);
  }

  list<point> CH=CONVEX_HULL(L);
  std::cout << "convex hull: ";
	
  while (!CH.empty()) {
    point p=CH.pop();	
    std::cout << "(" << p.xcoord() << "," << p.ycoord() << ")";
  }
  std::cout << endl;
	
  return 0;
}

Braided Lines (Verzopfte Geraden)

l1: y=9833x/9454   l2:y=9366x/9005.

Both lines path through the origin and the slope of l1 is slightly larger than the slope of l2.
The following program runs multiples of 0.001 between 0 and 1 and computes the corresponding y-values and, if the outcome of the comparison is different than in the previous iteration, prints x together with the current outcome.

The program outputs that l1 and l2 intersect several times, which is obviously not correct.

#include <iostream.h>

int main()
{
  cout.precision(12);
  float delta=0.001f;
  int last_comp=-1;
  float a=9833,b=9454,c=9366,d=9005;
  
  float x;
  for (x=0;x<0.1;x=x+delta) {
    float y1=a*x/b;
    float y2=c*x/d;
	
    int comp;
    if (y1<y2) comp=-1;
    else if (y1==y2) comp=0;
    else comp=1;
	
    if (comp!=last_comp) {
      cout << endl << x << ": ";
      if (comp==-1) cout << "l1 is below l2";
      if (comp==0) cout << "l1 intersects l2";
      else cout << "l1 is above l2";
    }
	
    last_comp=comp;
  }
		
  cout << endl << endl;
  return 0;
}

See also:

Data Types for 2-Dimensional Geometry

Linear Lists

Number Types

Convex Hulls

Geometry

Geometry Algorithms