#include <iostream>
#include <vector>

#include <boost/multiprecision/gmp.hpp>

//#include <boost/multiprecision/cpp_int.hpp>
//#include <boost/multiprecision/number.hpp>
#include <boost/numeric/conversion/converter.hpp>

typedef boost::multiprecision::mpq_rational GeometryCoord;
typedef boost::numeric::conversion_traits<double , GeometryCoord> MyFloat_to_double_Traits;

namespace boost {
  namespace numeric {
template<>
struct raw_converter<MyFloat_to_double_Traits>
{
  static double low_level_convert ( GeometryCoord const&  s ) { return( 0.0 ); } // Must be real conversion
};
}
}

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/io/wkt/wkt.hpp>

#include <boost/foreach.hpp>

//typedef boost::multiprecision::number<boost::multiprecision::gmp_rational, boost::multiprecision::et_off> GeometryCoord;
//typedef boost::multiprecision::number<boost::multiprecision::gmp_rational> GeometryCoord;
//typedef boost::multiprecision::mpq_rational GeometryCoord;
//typedef boost::multiprecision::number<boost::multiprecision::cpp_rational, boost::multiprecision::et_off> GeometryCoord;
typedef boost::geometry::model::d2::point_xy<GeometryCoord> GeometryPoint;
typedef boost::geometry::model::segment<GeometryPoint> GeometrySegment;
typedef boost::geometry::model::polygon<GeometryPoint> GeometryPolygon;

int main()
{
    /*std::vector<GeometryPoint> points;

    GeometrySegment Segment1(GeometryPoint(0, 3), GeometryPoint(3, 2));
    GeometrySegment Segment2(GeometryPoint(2, 0), GeometryPoint(0, 5));

    boost::geometry::intersection(Segment1, Segment2, points);*/
/*
    for (std::vector<point>::const_iterator i = points.begin(); i != points.end(); ++i) {
        point p = *i;

	std::cout << "x = " << p.x() << " y = " << p.y() << std::endl;
    }

    std::cout << std::endl;*/

    GeometryCoord c(3, 4);

    std::cout << " c = " << c << std::endl;

    GeometryPolygon green, blue;

    /*boost::geometry::read_wkt(
        "POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))", green);   

    boost::geometry::read_wkt(
    "POLYGON((10 10, 10 20, 20 20, 20 10, 10 10))", blue);*/




    /*boost::geometry::read_wkt(
        "POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))", green);   

    boost::geometry::read_wkt(
    "POLYGON((0 0, 0 20, 20 20, 20 0, 0 0))", blue);*/

    /*boost::geometry::read_wkt(
        "POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))", green);   
	
    boost::geometry::read_wkt(
    "POLYGON((1 1, 1 2, 2 2, 2 1, 1 1))", blue);*/

    green.outer().push_back(GeometryPoint(GeometryCoord(0), GeometryCoord(0)));
    green.outer().push_back(GeometryPoint(GeometryCoord(0, 1), GeometryCoord(10, 1)));
    green.outer().push_back(GeometryPoint(GeometryCoord(101, 10), GeometryCoord(101, 10)));
    green.outer().push_back(GeometryPoint(GeometryCoord(10, 1), GeometryCoord(0, 1)));
    green.outer().push_back(GeometryPoint(GeometryCoord(0), GeometryCoord(0)));


    blue.outer().push_back(GeometryPoint(GeometryCoord(10, 1), GeometryCoord(10, 1)));
    blue.outer().push_back(GeometryPoint(GeometryCoord(10, 1), GeometryCoord(20, 1)));
    blue.outer().push_back(GeometryPoint(GeometryCoord(20, 1), GeometryCoord(20, 1)));
    blue.outer().push_back(GeometryPoint(GeometryCoord(20, 1), GeometryCoord(10, 1)));
    blue.outer().push_back(GeometryPoint(GeometryCoord(10, 1), GeometryCoord(10, 1)));


    std::vector<GeometryPolygon> output;
    //boost::geometry::union_(green, blue, output);

    boost::geometry::intersection(green, blue, output);

    //boost::geometry::correct(green);

    std::cout << "size = " << output.size() << std::endl;

    //std::cout << "Intersects = " << boost::geometry::intersects(green, blue) << std::endl;

    BOOST_FOREACH(GeometryPolygon const& p, output) {
      //std::cout << boost::geometry::wkt<polygon>(p) << std::endl;

      const boost::geometry::model::ring<GeometryPoint> &ring = p.outer();

      for(boost::geometry::model::ring<GeometryPoint>::const_iterator i = ring.begin(); i != ring.end(); ++i ) {
        GeometryPoint Point = *i;

	std::cout << "x = " << Point.x() << " y = " << Point.y() << std::endl;
      }

      //std::cout << "Point inside: " << boost::geometry::within( point(0, 0), p ) << std::endl;
    }

    return 0;
}