Verifying the correctness of side_robust

Readme.txt

The program below is meant to verify the correctness of the orientation predicate (side strategy) implemented in
boost/geometry/extensions/triangulation/strategies/cartesian/side_robust.hpp .

It depends on Boost, CGAL, MPFR and GMP. Because it uses a Boost.Geometry extension, it requires the develop tree of
Boost.Geometry (at commit 10ecf6565b64800f6935fea81216564b2bbe1598 at the time of writing this).

I compiled it with the command:

g++ -DNDEBUG -lmpfr -lgmp -lCGAL -I/path/to/boost/develop/tree/include -O3 -march=haswell main.cpp

The output was

continuous random coordinates.
Robust consistent? 1
Approx consistent? 1
Robust correct? 1
Approx correct? 1

continuous random coordinates, continuous random exponent.
Robust consistent? 1
Approx consistent? 0
Robust correct? 1
Approx correct? 0

grid coordinates, random perturbation.
Robust consistent? 1
Approx consistent? 0
Robust correct? 1
Approx correct? 0

main.cpp

#include <iostream>
#include <random>

#include <boost/geometry.hpp>

#include <boost/geometry/extensions/triangulation/strategies/cartesian/side_robust.hpp>

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>

typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef Kernel::Point_2 Point_2;

int sign(double a) { return a > 0 ? 1 : (a < 0 ? -1 : 0); }

template <typename Strategy, typename Point>
bool consistent(Point const &p1, Point const &p2, Point const &p3) {
  return sign(Strategy::apply(p1, p2, p3)) ==
             sign(Strategy::apply(p3, p1, p2)) &&
         sign(Strategy::apply(p1, p2, p3)) ==
             sign(Strategy::apply(p2, p3, p1)) &&
         sign(Strategy::apply(p1, p2, p3)) ==
             -sign(Strategy::apply(p2, p1, p3)) &&
         sign(Strategy::apply(p1, p2, p3)) ==
             -sign(Strategy::apply(p1, p3, p2)) &&
         sign(Strategy::apply(p1, p2, p3)) ==
             -sign(Strategy::apply(p3, p2, p1));
}

int main() {
  namespace bg = boost::geometry;
  using p = bg::model::d2::point_xy<double>;
  using robust = bg::strategy::side::side_robust<double>;
  using approx = bg::strategy::side::side_by_triangle<double>;

  const int tests = 1'000'000;
  std::random_device rd;
  std::mt19937 gen(rd());
  bool consistent_approx = true, consistent_robust = true;
  bool correct_approx = true, correct_robust = true;

  // continuous random coordinates
  for (int i = 0; i < tests; ++i) {
    std::uniform_real_distribution<> dis(-40.0, 40.0);
    p p1(dis(gen), dis(gen));
    p p2(dis(gen), dis(gen));
    p p3(dis(gen), dis(gen));

    Point_2 cp1(p1.x(), p1.y()), cp2(p2.x(), p2.y()), cp3(p3.x(), p3.y());
    auto det = CGAL::orientation(cp1, cp2, cp3);
    int csign = det > 0 ? 1 : (det < 0 ? -1 : 0);
    if (!consistent<robust>(p1, p2, p3))
      consistent_robust = false;
    if (!consistent<approx>(p1, p2, p3))
      consistent_approx = false;
    if (csign != sign(robust::apply(p1, p2, p3)))
      correct_robust = false;
    if (csign != sign(approx::apply(p1, p2, p3)))
      correct_approx = false;
  }
  std::cout << "continuous random coordinates.\n";
  std::cout << "Robust consistent? " << consistent_robust << "\n";
  std::cout << "Approx consistent? " << consistent_approx << "\n";
  std::cout << "Robust correct? " << correct_robust << "\n";
  std::cout << "Approx correct? " << correct_approx << "\n\n";
  // continuous random coordinates, continuous random exponent
  consistent_approx = consistent_robust = correct_approx = correct_robust =
      true;
  for (int i = 0; i < tests; ++i) {
    std::uniform_real_distribution<> dis(-40.0, 40.0);
    p p1(dis(gen) * std::pow(10.0, dis(gen)),
         dis(gen) * std::pow(10.0, dis(gen)));
    p p2(dis(gen) * std::pow(10.0, dis(gen)),
         dis(gen) * std::pow(10.0, dis(gen)));
    p p3(dis(gen) * std::pow(10.0, dis(gen)),
         dis(gen) * std::pow(10.0, dis(gen)));

    Point_2 cp1(p1.x(), p1.y()), cp2(p2.x(), p2.y()), cp3(p3.x(), p3.y());
    auto det = CGAL::orientation(cp1, cp2, cp3);
    int csign = det > 0 ? 1 : (det < 0 ? -1 : 0);
    if (!consistent<robust>(p1, p2, p3))
      consistent_robust = false;
    if (!consistent<approx>(p1, p2, p3))
      consistent_approx = false;
    if (csign != sign(robust::apply(p1, p2, p3)))
      correct_robust = false;
    if (csign != sign(approx::apply(p1, p2, p3)))
      correct_approx = false;
  }
  std::cout << "continuous random coordinates, continuous random exponent.\n";
  std::cout << "Robust consistent? " << consistent_robust << "\n";
  std::cout << "Approx consistent? " << consistent_approx << "\n";
  std::cout << "Robust correct? " << correct_robust << "\n";
  std::cout << "Approx correct? " << correct_approx << "\n\n";

  // grid coordinates, random perturbation
  consistent_approx = consistent_robust = correct_approx = correct_robust =
      true;
  for (int i = 0; i < tests; ++i) {
    std::uniform_int_distribution<> dis(-40, 40);
    p p1(dis(gen), dis(gen));
    p p2(dis(gen), dis(gen));
    p p3(dis(gen), dis(gen));

    namespace trans = bg::strategy::transform;
    if (dis(gen) > 0) {
      std::uniform_real_distribution<> disa(0.0, 360.0);
      trans::rotate_transformer<bg::degree, double, 2, 2> rotate(disa(gen));
      p rp1, rp2, rp3;
      bg::transform(p1, rp1, rotate);
      bg::transform(p2, rp2, rotate);
      bg::transform(p3, rp3, rotate);
      p1 = rp1;
      p2 = rp2;
      p3 = rp3;
    }
    if (dis(gen) > 0) {
      std::uniform_real_distribution<> diss(-40.0, 40.0);
      trans::scale_transformer<double, 2, 2> scale(diss(gen));
      p sp1, sp2, sp3;
      bg::transform(p1, sp1, scale);
      bg::transform(p2, sp2, scale);
      bg::transform(p3, sp3, scale);
      p1 = sp1;
      p2 = sp2;
      p3 = sp3;
    }

    Point_2 cp1(p1.x(), p1.y()), cp2(p2.x(), p2.y()), cp3(p3.x(), p3.y());
    auto det = CGAL::orientation(cp1, cp2, cp3);
    int csign = det > 0 ? 1 : (det < 0 ? -1 : 0);
    if (!consistent<robust>(p1, p2, p3))
      consistent_robust = false;
    if (!consistent<approx>(p1, p2, p3))
      consistent_approx = false;
    if (csign != sign(robust::apply(p1, p2, p3)))
      correct_robust = false;
    if (csign != sign(approx::apply(p1, p2, p3)))
      correct_approx = false;
  }
  std::cout << "grid coordinates, random perturbation.\n";
  std::cout << "Robust consistent? " << consistent_robust << "\n";
  std::cout << "Approx consistent? " << consistent_approx << "\n";
  std::cout << "Robust correct? " << correct_robust << "\n";
  std::cout << "Approx correct? " << correct_approx << "\n\n";

  return 0;
}