Incorrect subtract result using Boost Polygon

I have the following two input polygons for which I want to calculate a subtracted polygon:

               * (0, 8)
              / \
             /   \
            /     \
   (-3, 0) *-------* (3, 0)

      (-1, 2) *-----* (1, 2)
              |     |
      (-1, 1) *-----* (1, 1)

Thus, I want to calculate A - B, which should result in a triangle with a square cutout. Calculating this using Boost Polygon results in an incorrect partial triangle with a cutout. It is hard to draw; the missing part of the result triangle is represented by the triangle (3, 0) => (0, 8) => (1, 2). I am using the following code to calculate the subtraction:

#include <boost/polygon/polygon.hpp>

namespace bp = boost::polygon;

int main()
{
  using Polygon = bp::polygon_data<int>;
  using Point = bp::point_data<int>;
  using PolygonSet = bp::polygon_set_data<int>;
  using SimplePolygons = std::vector<bp::polygon_data<int>>;

  using namespace boost::polygon::operators;

  Polygon A;
  {
    std::vector<Point> points{{-3, 0}, {3, 0}, {0, 8}};
    bp::set_points(A, points.begin(), points.end());
  }

  Polygon B;
  {
    std::vector<Point> points{{-1, 1}, {1, 1}, {1, 2}, {-1, 2}};
    bp::set_points(B, points.begin(), points.end());
  }

  PolygonSet result{A - B};

  SimplePolygons simplePolygons;
  result.get<SimplePolygons>(simplePolygons);
  for (const auto& polygon : simplePolygons)
  {
    for (const Point& p : polygon)
    {
      std::cout << '(' << std::to_string(p.x()) << ", " << std::to_string(p.y()) << ")\n";
    }
  }

  return 0;
}

This prints the following subsequent points making up the cutout triangle:

(3, 0)
(1, 2)
(1, 1)
(-1, 1)
(-1, 2)
(1, 2)
(0, 8)
(-3, 0)
(3, 0)

So, the edges (1, 2) => (3, 0) and (3, 0) => (0, 8) are missing from the result. The upper right part of the input triangle is missing from the result.

Correct output might look as follows:

(3, 0)
(1, 2)
(1, 1)
(-1, 1)
(-1, 2)
(1, 2)
(3, 0)
(0, 8)
(-3, 0)
(3, 0)

Is this a bug in Boost Polygon, am I somehow using the library incorrectly, or is there something else I am missing?

Some additional information:

I am using GCC 7.2.0, Boost 1.60.0. Boost Polygon has not been updated since, so upgrading Boost most likely will not help.
Parameterizing the point type and all other geometric types using double instead of int does not fix the issue.
Calculating cutouts using axis-aligned rectangles for example works without problems.
For my application I want to use Boost Polygon instead of Boost Geometry because it provides polygon keyhole fracturing support.

Solution

Answering my own question...

Boost Polygon was written with integer data types in mind. From the documentation:

In general a data type should define std::numeric_limits and be integer-like. Floating point coordinate types are not supported by all the algorithms and generally not suitable for use with the library at present (http://www.boost.org/doc/libs/1_60_0/libs/polygon/doc/gtl_coordinate_concept.htm)

I suspected this is some kind of precision issue I do not understand fully. Indeed, scaling all input coordinates by a factor of 1000 for example does result in a correct polygon:

(3000, 0)
(1000, 5333)
(1000, 2000)
(1000, 1000)
(-1000, 1000)
(-1000, 2000)
(1000, 2000)
(1000, 5333)
(0, 8000)
(-3000, 0)
(3000, 0)

So for the original input, it seems the keyhole fracturing algorithm is intent on adding a new vertex on the edge (3, 0) -> (0, 8) from which to enter the 'keyhole polygon'. The best possible vertex on integer grid position it can create for this is at (0, 8). So the result represents an approximation.

Indeed, providing the algorithm with similar input, for which a good candidate vertex exists on the triangle edge does result in the correct output. One such input triangle would be for example (-4, 0) - (4, 0) - (0, 8).

I see this as a limitation of the keyhole fracturing algorithm.