diff --git a/test/formulas/Jamfile.v2 b/test/formulas/Jamfile.v2
index dfd4b1878..18587bc24 100644
--- a/test/formulas/Jamfile.v2
+++ b/test/formulas/Jamfile.v2
@@ -12,5 +12,5 @@ test-suite boost-geometry-formulas
     :
     [ run inverse.cpp                        : : : : formulas_inverse ]
     [ run direct.cpp                         : : : : formulas_direct ]
-    #[ run intersection.cpp                   : : : : formulas_intersection ]
+    [ run intersection.cpp                   : : : : formulas_intersection ]
     ;
diff --git a/test/formulas/intersection.cpp b/test/formulas/intersection.cpp
new file mode 100644
index 000000000..f0b1b49d7
--- /dev/null
+++ b/test/formulas/intersection.cpp
@@ -0,0 +1,87 @@
+// Boost.Geometry
+// Unit Test
+
+// Copyright (c) 2016 Oracle and/or its affiliates.
+
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#include "test_formula.hpp"
+#include "intersection_cases.hpp"
+
+#include <boost/geometry/formulas/andoyer_inverse.hpp>
+#include <boost/geometry/formulas/gnomonic_intersection.hpp>
+#include <boost/geometry/formulas/sjoberg_intersection.hpp>
+#include <boost/geometry/formulas/thomas_direct.hpp>
+#include <boost/geometry/formulas/thomas_inverse.hpp>
+#include <boost/geometry/formulas/vincenty_direct.hpp>
+#include <boost/geometry/formulas/vincenty_inverse.hpp>
+
+void check_inverse(expected_result const& result, expected_result const& expected, expected_result const& reference, double reference_error)
+{
+    check_one(result.lon, expected.lon, reference.lon, reference_error);
+    check_one(result.lat, expected.lat, reference.lat, reference_error);
+}
+
+void test_all(expected_results const& results)
+{
+    double const d2r = bg::math::d2r<double>();
+    double const r2d = bg::math::r2d<double>();
+
+    double lona1r = results.p1.lon * d2r;
+    double lata1r = results.p1.lat * d2r;
+    double lona2r = results.p2.lon * d2r;
+    double lata2r = results.p2.lat * d2r;
+    double lonb1r = results.q1.lon * d2r;
+    double latb1r = results.q1.lat * d2r;
+    double lonb2r = results.q2.lon * d2r;
+    double latb2r = results.q2.lat * d2r;
+
+    expected_result result;
+
+    // WGS84
+    bg::srs::spheroid<double> spheroid(6378137.0, 6356752.3142451793);
+
+    bg::formula::gnomonic_intersection<double, bg::formula::vincenty_inverse, bg::formula::vincenty_direct>
+        ::apply(lona1r, lata1r, lona2r, lata2r, lonb1r, latb1r, lonb2r, latb2r, result.lon, result.lat, spheroid);
+    result.lon *= r2d;
+    result.lat *= r2d;
+    check_inverse(result, results.gnomonic_vincenty, results.gnomonic_karney, 0.00000001);
+
+    bg::formula::gnomonic_intersection<double, bg::formula::thomas_inverse, bg::formula::thomas_direct>
+        ::apply(lona1r, lata1r, lona2r, lata2r, lonb1r, latb1r, lonb2r, latb2r, result.lon, result.lat, spheroid);
+    result.lon *= r2d;
+    result.lat *= r2d;
+    check_inverse(result, results.gnomonic_thomas, results.gnomonic_karney, 0.0000001);
+
+    bg::formula::sjoberg_intersection<double, bg::formula::vincenty_inverse, 4>
+        ::apply(lona1r, lata1r, lona2r, lata2r, lonb1r, latb1r, lonb2r, latb2r, result.lon, result.lat, spheroid);
+    result.lon *= r2d;
+    result.lat *= r2d;
+    check_inverse(result, results.sjoberg_vincenty, results.sjoberg_karney, 0.00000001);
+
+    bg::formula::sjoberg_intersection<double, bg::formula::thomas_inverse, 2>
+        ::apply(lona1r, lata1r, lona2r, lata2r, lonb1r, latb1r, lonb2r, latb2r, result.lon, result.lat, spheroid);
+    result.lon *= r2d;
+    result.lat *= r2d;
+    check_inverse(result, results.sjoberg_thomas, results.sjoberg_karney, 0.0000001);
+
+    bg::formula::sjoberg_intersection<double, bg::formula::andoyer_inverse, 1>
+        ::apply(lona1r, lata1r, lona2r, lata2r, lonb1r, latb1r, lonb2r, latb2r, result.lon, result.lat, spheroid);
+    result.lon *= r2d;
+    result.lat *= r2d;
+    check_inverse(result, results.sjoberg_andoyer, results.sjoberg_karney, 0.0001);
+}
+
+int test_main(int, char*[])
+{
+    for (size_t i = 0; i < expected_size; ++i)
+    {
+        test_all(expected[i]);
+    }
+
+    return 0;
+}
diff --git a/test/formulas/intersection_cases.hpp b/test/formulas/intersection_cases.hpp
new file mode 100644
index 000000000..c0c40d3a4
--- /dev/null
+++ b/test/formulas/intersection_cases.hpp
@@ -0,0 +1,129 @@
+// Boost.Geometry
+// Unit Test
+
+// Copyright (c) 2016 Oracle and/or its affiliates.
+
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_TEST_INTERSECTION_CASES_HPP
+#define BOOST_GEOMETRY_TEST_INTERSECTION_CASES_HPP
+
+struct coordinates
+{
+    double lon;
+    double lat;
+};
+
+struct expected_result
+{
+    double lon;
+    double lat;
+};
+
+struct expected_results
+{
+    coordinates p1;
+    coordinates p2;
+    coordinates q1;
+    coordinates q2;
+    expected_result gnomonic_karney;
+    expected_result gnomonic_vincenty;
+    expected_result gnomonic_thomas;
+    expected_result sjoberg_karney;
+    expected_result sjoberg_vincenty;
+    expected_result sjoberg_thomas;
+    expected_result sjoberg_andoyer;
+};
+
+expected_results expected[] =
+{
+    {
+        { -1, -1 },{ 1, 1 },
+        { -1, 1 },{ 1, -1 },
+        { 0.0000000000000000, 0.0000000000000002 },
+        { 0.0000000000000000, 0.0000000000000000 },
+        { 0.0000000000000000, 0.0000000000000000 },
+        { 0.0000000000633173, 0.0000000000000003 },
+        { 0.0000000000626537, -0.0000000000000000 },
+        { -0.0000000187861002, -0.0000000000000000 },
+        { -0.0000055778585615, -0.0000000000000000 }
+    },{
+        { 1, 1 },{ -1, -1 },
+        { -1, 1 },{ 1, -1 },
+        { 0.0000000000000000, 0.0000000000000002 },
+        { 0.0000000000000000, 0.0000000000000000 },
+        { 0.0000000000000000, 0.0000000000000000 },
+        { 0.0000000000000000, 0.0000000000633274 },
+        { 0.0000000000000000, 0.0000000000626632 },
+        { -0.0000000000000006, -0.0000000187889745 },
+        { -0.0000000000000001, -0.0000055787431353 }
+    },{
+        { -1, -1 },{ 1, 1 },
+        { 1, -1 },{ -1, 1 },
+        { 0.0000000000000000, 0.0000000000000002 },
+        { 0.0000000000000000, 0.0000000000000000 },
+        { 0.0000000000000000, 0.0000000000000000 },
+        { -0.0000000000000000, -0.0000000000633268 },
+        { -0.0000000000000000, -0.0000000000626632 },
+        { -0.0000000000000004, 0.0000000187889746 },
+        { 0.0000000000000001, 0.0000055787431353 }
+    },{
+        { 1, 1 },{ -1, -1 },
+        { 1, -1 },{ -1, 1 },
+        { 0.0000000000000000, 0.0000000000000002 },
+        { 0.0000000000000000, 0.0000000000000000 },
+        { 0.0000000000000000, 0.0000000000000000 },
+        { -0.0000000000633173, 0.0000000000000003 },
+        { -0.0000000000626537, -0.0000000000000000 },
+        { 0.0000000187860994, 0.0000000000000001 },
+        { 0.0000055778585615, -0.0000000000000000 }
+    },{
+        { 0, 0 },{ 1, 1 },
+        { 0, 1 },{ 1, 0 },
+        { 0.5000000000000000, 0.5000573755188484 },
+        { 0.5000000000000000, 0.5000573755188470 },
+        { 0.5000000000000000, 0.5000573755109839 },
+        { 0.5000000000316606, 0.5000573755188390 },
+        { 0.5000000000313266, 0.5000573755188389 },
+        { 0.4999999906069524, 0.5000573755152582 },
+        { 0.4999972102164753, 0.5000573755151276 }
+    },{
+        { 1, 1 },{ 0, 0 },
+        { 0, 1 },{ 1, 0 },
+        { 0.5000000000000000, 0.5000573755188484 },
+        { 0.5000000000000000, 0.5000573755188470 },
+        { 0.5000000000000000, 0.5000573755109839 },
+        { 0.5000000000000000, 0.5000573755505008 },
+        { 0.5000000000000000, 0.5000573755501669 },
+        { 0.4999999999999996, 0.5000573661218464 },
+        { 0.4999999999999999, 0.5000545856093679 }
+    },{
+        { 0, 0 },{ 1, 1 },
+        { 1, 0 },{ 0, 1 },
+        { 0.5000000000000000, 0.5000573755188484 },
+        { 0.5000000000000000, 0.5000573755188470 },
+        { 0.5000000000000000, 0.5000573755109839 },
+        { 0.5000000000000001, 0.5000573754871773 },
+        { 0.4999999999999999, 0.5000573754875109 },
+        { 0.4999999999999999, 0.5000573849086647 },
+        { 0.5000000000000000, 0.5000601654208935 }
+    },{
+        { 1, 1 },{ 0, 0 },
+        { 1, 0 },{ 0, 1 },
+        { 0.5000000000000000, 0.5000573755188484 },
+        { 0.5000000000000000, 0.5000573755188470 },
+        { 0.5000000000000000, 0.5000573755109839 },
+        { 0.4999999999683394, 0.5000573755188390 },
+        { 0.4999999999686731, 0.5000573755188389 },
+        { 0.5000000093930521, 0.5000573755152582 },
+        { 0.5000027897835244, 0.5000573755151276 }
+    }
+};
+
+size_t expected_size = sizeof(expected) / sizeof(expected_results);
+
+#endif // BOOST_GEOMETRY_TEST_INTERSECTION_CASES_HPP
diff --git a/test/formulas/test_formula.hpp b/test/formulas/test_formula.hpp
index f16dbd79f..92f690bd3 100644
--- a/test/formulas/test_formula.hpp
+++ b/test/formulas/test_formula.hpp
@@ -34,13 +34,17 @@ void check_one(double result, double expected, double reference, double referenc
     }
 
     double res_max = (std::max)(bg::math::abs(result), bg::math::abs(expected));
-    if (res_max > 10 * std::numeric_limits<double>::epsilon())
+    if (res_max > 100 * std::numeric_limits<double>::epsilon())
     {
         BOOST_CHECK_CLOSE(result, expected, 0.001);
     }
+    else if (res_max > 10 * std::numeric_limits<double>::epsilon())
+    {
+        BOOST_CHECK_CLOSE(result, expected, 0.1);
+    }
     else if (res_max > std::numeric_limits<double>::epsilon())
     {
-        BOOST_CHECK_CLOSE(result, expected, 0.01);
+        BOOST_CHECK_CLOSE(result, expected, 10);
     }
 
     // NOTE: in some cases it probably will be necessary to normalize