c++ - How to use Boost::Geometry _union with integers

Question

I'm trying to use Boost::Geometry _union with integers, for both performance and numeric accuracy. For that, I multiply the coordinates of my input with 10,000. Thus creating coordinates with up to 9 digits. I figured that since I use 64-bit integers, that should work well.

Unfortunately, when I run the code, I get strange results (the output polygon includes a point that is far from any polygon in the input). Investigating the code of Boost::Geometry brought me to the conclusion that the origin is a wrap-around problem in the file cart_intersect.hpp:

set<1>(point, get<0, 1>(segment) + boost::numeric_cast
            <
                CoordinateType
            >(numerator * dy_promoted / denominator));

When all three (numerator,dy_promoted and denominator) are huge, the result of the multiplication takes more than 64-bit, and therefore the whole result is incorrect.

Is the use of Boost::Geometry with such big integers valid? What is the correct way to use Boost::Geometry with integers while maintaining correctness and precision?

---

Edit @sehe, thanks for your response. Here is an SSCCE, compiled with VS2013 and Boost 1.63

#include <boost/geometry/geometry.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/multi/geometries/multi_polygon.hpp> 
#include <iostream>
#include <string>
#include <vector>

using namespace std;
namespace bg = boost::geometry;

typedef bg::model::d2::point_xy<long long, bg::cs::cartesian> TPoint;
typedef bg::model::ring<TPoint> TRing;
typedef bg::model::polygon<TPoint> TPolygon;
typedef bg::model::multi_polygon<TPolygon> TMultiPolygon;

void PrintRing(const TRing& rng)
{
  for (const auto& ver : rng)
  {
    cout << "(" << ver.x() << "," << ver.y() << "),";
  }
}

void PrintPolygon(const TPolygon& pol)
{
  cout << "Outer: ";
  PrintRing(pol.outer());
  cout << endl;

  for (const auto& rng : pol.inners())
  {
    cout << "Inner: ";
    PrintRing(rng);
    cout << endl;
  }
}

void PrintMultiPolygon(const string name, const TMultiPolygon& mp)
{
  cout << "Multi-Polygon " << name << " : " << endl;
  for (const auto& pol : mp)
  {
    PrintPolygon(pol);
  }
  cout << endl;
}

int main()
{
  cout << "BOOST_LIB_VERSION: " << BOOST_LIB_VERSION << endl;
  const vector<TPoint> verticesA{ { -405129, 2010409 }, { 3370580, 2010409 }, { 3370580, 1997709 }, { -405129, 1997709 }, { -405129, 2010409 } };
  const TRing rngA(verticesA.cbegin(), verticesA.cend());
  TPolygon polA;
  polA.outer() = rngA;
  TMultiPolygon mpA;
  mpA.push_back(polA);

  const vector<TPoint> verticesB{ { 3364230, -895349 }, { 3364230, 2004060 }, { 3376930, 2004059 }, { 3376930, -895350 }, { 3364230, -895349 } };
  const TRing rngB(verticesB.cbegin(), verticesB.cend());
  TPolygon polB;
  polB.outer() = rngB;
  TMultiPolygon mpB;
  mpB.push_back(polB);

  TMultiPolygon output;

  bg::union_(mpA, mpB, output);

  PrintMultiPolygon("A", mpA);
  PrintMultiPolygon("B", mpB);
  PrintMultiPolygon("output", output);
}

---

The output of the program is:

BOOST_LIB_VERSION: 1_63

Multi-Polygon A :
Outer: (-405129,2010409),(3370580,2010409),(3370580,1997709),(-405129,1997709),(-405129,2010409),

Multi-Polygon B :
Outer: (3364230,-895349),(3364230,2004060),(3376930,2004059),(3376930,-895350),(3364230,-895349),

Multi-Polygon output :
Outer: (3370580,2004060),(3376930,2004059),(3376930,-895350),(3364230,-895349),(3364230,-1372382),(-405129,1997709),(-405129,2010409),(3370580,2010409),(3370580,2004060),

Note the coordinates in bold, the Y value is farther than any Y coordinate in the input.

Answer 1

Indeed, running with -fsanitize=undefined prints

 /home/sehe/custom/boost/boost/geometry/strategies/cartesian/intersection.hpp:190:18: runtime error: signed integer overflow: 10923345128122 * 2899409 cannot be represented in type 'long long int'

Instead you could use a vendor-defined 128 bit extension, or use Boost's:

namespace mp = boost::multiprecision;

typedef mp::checked_int128_t T;
typedef bg::model::d2::point_xy<T, bg::cs::cartesian> TPoint;

In fact, you can use a arbitrary-precision integer:

typedef mp::checked_cpp_int T;

Note If you want to use unchecked arithmetic with cpp_int you will have to make sure that expression-templates are disabled for boost

typedef mp::number<mp::backends::cpp_int_backend<0, 0, mp::signed_magnitude, mp::unchecked>, mp::et_off> T;

See e.g. Why does using boost::multiprecision::cpp_int affect tail call optimization here, How to use sqrt and ceil with Boost::multiprecision?, etc.

With all the above the output becomes:

Live On Wandbox

BOOST_LIB_VERSION: 1_64
Multi-Polygon A : 
Outer: (-405129,2010409),(3370580,2010409),(3370580,1997709),(-405129,1997709),(-405129,2010409),

Multi-Polygon B : 
Outer: (3364230,-895349),(3364230,2004060),(3376930,2004059),(3376930,-895350),(3364230,-895349),

Multi-Polygon output : 
Outer: (3370580,2004060),(3376930,2004059),(3376930,-895350),(3364230,-895349),(3364230,1997709),(-405129,1997709),(-405129,2010409),(3370580,2010409),(3370580,2004060),

See more: Boost Geometry and exact point types