#include <iomanip>
 #include <iostream>
 #include <vector>
 #include <boost/geometry.hpp>
 #include <boost/geometry/geometries/point_xy.hpp>
 #include <boost/geometry/geometries/polygon.hpp>
 #include <boost/geometry/geometries/register/point.hpp>
 #include <boost/geometry/geometries/register/ring.hpp>
 #include <boost/geometry/io/wkt/wkt.hpp>
 #include <boost/geometry/multi/geometries/multi_polygon.hpp>

typedef boost::geometry::model::d2::point_xy<double> pt;
typedef boost::geometry::model::polygon<pt> polygon;
typedef boost::geometry::model::multi_polygon<polygon> multi_polygon;

int main()
{
    int num_orig=50;
    int num_rounds=100000;
    srand(1234);
    std::cout<<std::setprecision(16);

    std::vector<multi_polygon> poly_list;

    for(int i=0;i<num_orig;i++)
    {
        multi_polygon mp;
        polygon p;
        for(int j=0;j<3;j++)
        {
            double x=(double)rand()/RAND_MAX;
            double y=(double)rand()/RAND_MAX;
            p.outer().push_back(pt(x,y));
        }
        boost::geometry::correct(p);
        mp.push_back(p);
        boost::geometry::detail::overlay::has_self_intersections(mp);
        poly_list.push_back(mp);
    }


    for(int j=0;j<num_rounds;j++)
    {
        int a = rand() % poly_list.size();
        int b = rand() % poly_list.size();
        multi_polygon mp_i,mp_d,mp_e;
        boost::geometry::intersection(poly_list[a],poly_list[b],mp_i);
        boost::geometry::difference(poly_list[a],poly_list[b],mp_d);
        boost::geometry::difference(poly_list[b],poly_list[a],mp_e);
        try
        {
            boost::geometry::detail::overlay::has_self_intersections(mp_i);
        }
        catch(...)
        {
            std::cout<<"FAILED TO INTERSECT"<<std::endl;
            std::cout<<boost::geometry::wkt(poly_list[a])<<std::endl;
            std::cout<<boost::geometry::wkt(poly_list[b])<<std::endl;
            std::cout<<boost::geometry::wkt(mp_i)<<std::endl;
            break;
        }
        try
        {
            boost::geometry::detail::overlay::has_self_intersections(mp_d);
        }
        catch(...)
        {
            std::cout<<"FAILED TO SUBTRACT"<<std::endl;
            std::cout<<boost::geometry::wkt(poly_list[a])<<std::endl;
            std::cout<<boost::geometry::wkt(poly_list[b])<<std::endl;
            std::cout<<boost::geometry::wkt(mp_d)<<std::endl;
            break;
        }
        try
        {
            boost::geometry::detail::overlay::has_self_intersections(mp_e);
        }
        catch(...)
        {
            std::cout<<"FAILED TO SUBTRACT"<<std::endl;
            std::cout<<boost::geometry::wkt(poly_list[b])<<std::endl;
            std::cout<<boost::geometry::wkt(poly_list[a])<<std::endl;
            std::cout<<boost::geometry::wkt(mp_e)<<std::endl;
            break;
        }
        if(boost::geometry::area(mp_i)>0) poly_list.push_back(mp_i);
        if(boost::geometry::area(mp_d)>0) poly_list.push_back(mp_d);
        if(boost::geometry::area(mp_e)>0) poly_list.push_back(mp_e);
    }

    return 0;
}