[ceph.git] / ceph / src / boost / libs / geometry / include / boost / geometry / strategies / cartesian / side_of_intersection.hpp

// Boost.Geometry (aka GGL, Generic Geometry Library)

// Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands.

// This file was modified by Oracle on 2015.
// Modifications copyright (c) 2015, 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_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP
#define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP


#include <limits>

#include <boost/core/ignore_unused.hpp>
#include <boost/type_traits/is_integral.hpp>
#include <boost/type_traits/make_unsigned.hpp>

#include <boost/geometry/arithmetic/determinant.hpp>
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/assert.hpp>
#include <boost/geometry/core/coordinate_type.hpp>
#include <boost/geometry/algorithms/detail/assign_indexed_point.hpp>
#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
#include <boost/geometry/util/math.hpp>

#ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
#include <boost/math/common_factor_ct.hpp>
#include <boost/math/common_factor_rt.hpp>
#include <boost/multiprecision/cpp_int.hpp>
#endif

namespace boost { namespace geometry
{

namespace strategy { namespace side
{

namespace detail
{

// A tool for multiplication of integers avoiding overflow
// It's a temporary workaround until we can use Multiprecision
// The algorithm is based on Karatsuba algorithm
// see: http://en.wikipedia.org/wiki/Karatsuba_algorithm
template <typename T>
struct multiplicable_integral
{
    // Currently this tool can't be used with non-integral coordinate types.
    // Also side_of_intersection strategy sign_of_product() and sign_of_compare()
    // functions would have to be modified to properly support floating-point
    // types (comparisons and multiplication).
    BOOST_STATIC_ASSERT(boost::is_integral<T>::value);

    static const std::size_t bits = CHAR_BIT * sizeof(T);
    static const std::size_t half_bits = bits / 2;
    typedef typename boost::make_unsigned<T>::type unsigned_type;
    static const unsigned_type base = unsigned_type(1) << half_bits; // 2^half_bits

    int m_sign;
    unsigned_type m_ms;
    unsigned_type m_ls;

    multiplicable_integral(int sign, unsigned_type ms, unsigned_type ls)
        : m_sign(sign), m_ms(ms), m_ls(ls)
    {}

    explicit multiplicable_integral(T const& val)
    {
        unsigned_type val_u = val > 0 ?
                                unsigned_type(val)
                              : val == (std::numeric_limits<T>::min)() ?
                                  unsigned_type((std::numeric_limits<T>::max)()) + 1
                                : unsigned_type(-val);
        // MMLL -> S 00MM 00LL
        m_sign = math::sign(val);
        m_ms = val_u >> half_bits; // val_u / base
        m_ls = val_u - m_ms * base;
    }
    
    friend multiplicable_integral operator*(multiplicable_integral const& a,
                                            multiplicable_integral const& b)
    {
        // (S 00MM 00LL) * (S 00MM 00LL) -> (S Z2MM 00LL)
        unsigned_type z2 = a.m_ms * b.m_ms;
        unsigned_type z0 = a.m_ls * b.m_ls;
        unsigned_type z1 = (a.m_ms + a.m_ls) * (b.m_ms + b.m_ls) - z2 - z0;
        // z0 may be >= base so it must be normalized to allow comparison
        unsigned_type z0_ms = z0 >> half_bits; // z0 / base
        return multiplicable_integral(a.m_sign * b.m_sign,
                                      z2 * base + z1 + z0_ms,
                                      z0 - base * z0_ms);
    }

    friend bool operator<(multiplicable_integral const& a,
                          multiplicable_integral const& b)
    {
        if ( a.m_sign == b.m_sign )
        {
            bool u_less = a.m_ms < b.m_ms
                      || (a.m_ms == b.m_ms && a.m_ls < b.m_ls);
            return a.m_sign > 0 ? u_less : (! u_less);
        }
        else
        {
            return a.m_sign < b.m_sign;
        }
    }

    friend bool operator>(multiplicable_integral const& a,
                          multiplicable_integral const& b)
    {
        return b < a;
    }

    template <typename CmpVal>
    void check_value(CmpVal const& cmp_val) const
    {
        unsigned_type b = base; // a workaround for MinGW - undefined reference base
        CmpVal val = CmpVal(m_sign) * (CmpVal(m_ms) * CmpVal(b) + CmpVal(m_ls));
        BOOST_GEOMETRY_ASSERT(cmp_val == val);
    }
};

} // namespace detail

// Calculates the side of the intersection-point (if any) of
// of segment a//b w.r.t. segment c
// This is calculated without (re)calculating the IP itself again and fully
// based on integer mathematics; there are no divisions
// It can be used for either integer (rescaled) points, and also for FP
class side_of_intersection
{
private :
    template <typename T, typename U>
    static inline
    int sign_of_product(T const& a, U const& b)
    {
        return a == 0 || b == 0 ? 0
            : a > 0 && b > 0 ? 1
            : a < 0 && b < 0 ? 1
            : -1;
    }

    template <typename T>
    static inline
    int sign_of_compare(T const& a, T const& b, T const& c, T const& d)
    {
        // Both a*b and c*d are positive
        // We have to judge if a*b > c*d

        using side::detail::multiplicable_integral;
        multiplicable_integral<T> ab = multiplicable_integral<T>(a)
                                     * multiplicable_integral<T>(b);
        multiplicable_integral<T> cd = multiplicable_integral<T>(c)
                                     * multiplicable_integral<T>(d);
        
        int result = ab > cd ? 1
                   : ab < cd ? -1
                   : 0
                   ;

#ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
        using namespace boost::multiprecision;
        cpp_int const lab = cpp_int(a) * cpp_int(b);
        cpp_int const lcd = cpp_int(c) * cpp_int(d);

        ab.check_value(lab);
        cd.check_value(lcd);

        int result2 = lab > lcd ? 1
                    : lab < lcd ? -1
                    : 0
                    ;
        BOOST_GEOMETRY_ASSERT(result == result2);
#endif

        return result;
    }

    template <typename T>
    static inline
    int sign_of_addition_of_two_products(T const& a, T const& b, T const& c, T const& d)
    {
        // sign of a*b+c*d, 1 if positive, -1 if negative, else 0
        int const ab = sign_of_product(a, b);
        int const cd = sign_of_product(c, d);
        if (ab == 0)
        {
            return cd;
        }
        if (cd == 0)
        {
            return ab;
        }

        if (ab == cd)
        {
            // Both positive or both negative
            return ab;
        }

        // One is positive, one is negative, both are non zero
        // If ab is positive, we have to judge if a*b > -c*d (then 1 because sum is positive)
        // If ab is negative, we have to judge if c*d > -a*b (idem)
        return ab == 1
            ? sign_of_compare(a, b, -c, d)
            : sign_of_compare(c, d, -a, b);
    }


public :

    // Calculates the side of the intersection-point (if any) of
    // of segment a//b w.r.t. segment c
    // This is calculated without (re)calculating the IP itself again and fully
    // based on integer mathematics
    template <typename T, typename Segment, typename Point>
    static inline T side_value(Segment const& a, Segment const& b,
                Segment const& c, Point const& fallback_point)
    {
        // The first point of the three segments is reused several times
        T const ax = get<0, 0>(a);
        T const ay = get<0, 1>(a);
        T const bx = get<0, 0>(b);
        T const by = get<0, 1>(b);
        T const cx = get<0, 0>(c);
        T const cy = get<0, 1>(c);

        T const dx_a = get<1, 0>(a) - ax;
        T const dy_a = get<1, 1>(a) - ay;

        T const dx_b = get<1, 0>(b) - bx;
        T const dy_b = get<1, 1>(b) - by;

        T const dx_c = get<1, 0>(c) - cx;
        T const dy_c = get<1, 1>(c) - cy;

        // Cramer's rule: d (see cart_intersect.hpp)
        T const d = geometry::detail::determinant<T>
                    (
                        dx_a, dy_a,
                        dx_b, dy_b
                    );

        T const zero = T();
        if (d == zero)
        {
            // There is no IP of a//b, they are collinear or parallel
            // Assuming they intersect (this method should be called for
            // segments known to intersect), they are collinear and overlap.
            // They have one or two intersection points - we don't know and
            // have to rely on the fallback intersection point

            Point c1, c2;
            geometry::detail::assign_point_from_index<0>(c, c1);
            geometry::detail::assign_point_from_index<1>(c, c2);
            return side_by_triangle<>::apply(c1, c2, fallback_point);
        }

        // Cramer's rule: da (see cart_intersect.hpp)
        T const da = geometry::detail::determinant<T>
                    (
                        dx_b,    dy_b,
                        ax - bx, ay - by
                    );

        // IP is at (ax + (da/d) * dx_a, ay + (da/d) * dy_a)
        // Side of IP is w.r.t. c is: determinant(dx_c, dy_c, ipx-cx, ipy-cy)
        // We replace ipx by expression above and multiply each term by d

#ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
        T const result1 = geometry::detail::determinant<T>
                    (
                        dx_c * d,                   dy_c * d,
                        d * (ax - cx) + dx_a * da,  d * (ay - cy) + dy_a * da
                    );

        // Note: result / (d * d)
        // is identical to the side_value of side_by_triangle
        // Therefore, the sign is always the same as that result, and the
        // resulting side (left,right,collinear) is the same

        // The first row we divide again by d because of determinant multiply rule
        T const result2 = d * geometry::detail::determinant<T>
                    (
                        dx_c,                   dy_c,
                        d * (ax - cx) + dx_a * da,  d * (ay - cy) + dy_a * da
                    );
        // Write out:
        T const result3 = d * (dx_c * (d * (ay - cy) + dy_a * da)
                             - dy_c * (d * (ax - cx) + dx_a * da));
        // Write out in braces:
        T const result4 = d * (dx_c * d * (ay - cy) + dx_c * dy_a * da
                             - dy_c * d * (ax - cx) - dy_c * dx_a * da);
        // Write in terms of d * XX + da * YY
        T const result5 = d * (d * (dx_c * (ay - cy) - dy_c * (ax - cx))
                             + da * (dx_c * dy_a - dy_c * dx_a));

        boost::ignore_unused(result1, result2, result3, result4, result5);
        //return result;
#endif

        // We consider the results separately
        // (in the end we only have to return the side-value 1,0 or -1)

        // To avoid multiplications we judge the product (easy, avoids *d)
        // and the sign of p*q+r*s (more elaborate)
        T const result = sign_of_product
            (
                d,
                sign_of_addition_of_two_products
                    (
                        d, dx_c * (ay - cy) - dy_c * (ax - cx),
                        da, dx_c * dy_a - dy_c * dx_a
                    )
                );
        return result;


    }

    template <typename Segment, typename Point>
    static inline int apply(Segment const& a, Segment const& b,
            Segment const& c,
            Point const& fallback_point)
    {
        typedef typename geometry::coordinate_type<Segment>::type coordinate_type;
        coordinate_type const s = side_value<coordinate_type>(a, b, c, fallback_point);
        coordinate_type const zero = coordinate_type();
        return math::equals(s, zero) ? 0
            : s > zero ? 1
            : -1;
    }

};


}} // namespace strategy::side

}} // namespace boost::geometry


#endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP
Commit	Line	Data
7c673cae FG	1	// Boost.Geometry (aka GGL, Generic Geometry Library)
	2
	3	// Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands.
	4
	5	// This file was modified by Oracle on 2015.
	6	// Modifications copyright (c) 2015, Oracle and/or its affiliates.
	7
	8	// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
	9
	10	// Use, modification and distribution is subject to the Boost Software License,
	11	// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
	12	// http://www.boost.org/LICENSE_1_0.txt)
	13
	14	#ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP
	15	#define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP
	16
	17
	18	#include <limits>
	19
	20	#include <boost/core/ignore_unused.hpp>
	21	#include <boost/type_traits/is_integral.hpp>
	22	#include <boost/type_traits/make_unsigned.hpp>
	23
	24	#include <boost/geometry/arithmetic/determinant.hpp>
	25	#include <boost/geometry/core/access.hpp>
	26	#include <boost/geometry/core/assert.hpp>
	27	#include <boost/geometry/core/coordinate_type.hpp>
	28	#include <boost/geometry/algorithms/detail/assign_indexed_point.hpp>
	29	#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
	30	#include <boost/geometry/util/math.hpp>
	31
	32	#ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
	33	#include <boost/math/common_factor_ct.hpp>
	34	#include <boost/math/common_factor_rt.hpp>
	35	#include <boost/multiprecision/cpp_int.hpp>
	36	#endif
	37
	38	namespace boost { namespace geometry
	39	{
	40
	41	namespace strategy { namespace side
	42	{
	43
	44	namespace detail
	45	{
	46
	47	// A tool for multiplication of integers avoiding overflow
	48	// It's a temporary workaround until we can use Multiprecision
	49	// The algorithm is based on Karatsuba algorithm
	50	// see: http://en.wikipedia.org/wiki/Karatsuba_algorithm
	51	template <typename T>
	52	struct multiplicable_integral
	53	{
	54	// Currently this tool can't be used with non-integral coordinate types.
	55	// Also side_of_intersection strategy sign_of_product() and sign_of_compare()
	56	// functions would have to be modified to properly support floating-point
	57	// types (comparisons and multiplication).
	58	BOOST_STATIC_ASSERT(boost::is_integral<T>::value);
	59
	60	static const std::size_t bits = CHAR_BIT * sizeof(T);
	61	static const std::size_t half_bits = bits / 2;
	62	typedef typename boost::make_unsigned<T>::type unsigned_type;
	63	static const unsigned_type base = unsigned_type(1) << half_bits; // 2^half_bits
	64
65	int m_sign;
66	unsigned_type m_ms;
67	unsigned_type m_ls;
68
69	multiplicable_integral(int sign, unsigned_type ms, unsigned_type ls)
70	: m_sign(sign), m_ms(ms), m_ls(ls)
71	{}
72
73	explicit multiplicable_integral(T const& val)
74	{
75	unsigned_type val_u = val > 0 ?
76	unsigned_type(val)
77	: val == (std::numeric_limits<T>::min)() ?
78	unsigned_type((std::numeric_limits<T>::max)()) + 1
79	: unsigned_type(-val);
80	// MMLL -> S 00MM 00LL
81	m_sign = math::sign(val);
82	m_ms = val_u >> half_bits; // val_u / base
83	m_ls = val_u - m_ms * base;
84	}
85
86	friend multiplicable_integral operator*(multiplicable_integral const& a,
87	multiplicable_integral const& b)
88	{
89	// (S 00MM 00LL) * (S 00MM 00LL) -> (S Z2MM 00LL)
90	unsigned_type z2 = a.m_ms * b.m_ms;
91	unsigned_type z0 = a.m_ls * b.m_ls;
92	unsigned_type z1 = (a.m_ms + a.m_ls) * (b.m_ms + b.m_ls) - z2 - z0;
93	// z0 may be >= base so it must be normalized to allow comparison
94	unsigned_type z0_ms = z0 >> half_bits; // z0 / base
95	return multiplicable_integral(a.m_sign * b.m_sign,
96	z2 * base + z1 + z0_ms,
97	z0 - base * z0_ms);
98	}
99
100	friend bool operator<(multiplicable_integral const& a,
101	multiplicable_integral const& b)
102	{
103	if ( a.m_sign == b.m_sign )
104	{
105	bool u_less = a.m_ms < b.m_ms
106	\|\| (a.m_ms == b.m_ms && a.m_ls < b.m_ls);
107	return a.m_sign > 0 ? u_less : (! u_less);
108	}
109	else
110	{
111	return a.m_sign < b.m_sign;
112	}
113	}
114
115	friend bool operator>(multiplicable_integral const& a,
116	multiplicable_integral const& b)
117	{
118	return b < a;
119	}
120
121	template <typename CmpVal>
122	void check_value(CmpVal const& cmp_val) const
123	{
124	unsigned_type b = base; // a workaround for MinGW - undefined reference base
125	CmpVal val = CmpVal(m_sign) * (CmpVal(m_ms) * CmpVal(b) + CmpVal(m_ls));
126	BOOST_GEOMETRY_ASSERT(cmp_val == val);
127	}
128	};
129
130	} // namespace detail
131
132	// Calculates the side of the intersection-point (if any) of
133	// of segment a//b w.r.t. segment c
134	// This is calculated without (re)calculating the IP itself again and fully
135	// based on integer mathematics; there are no divisions
136	// It can be used for either integer (rescaled) points, and also for FP
137	class side_of_intersection
138	{
139	private :
140	template <typename T, typename U>
141	static inline
142	int sign_of_product(T const& a, U const& b)
143	{
144	return a == 0 \|\| b == 0 ? 0
145	: a > 0 && b > 0 ? 1
146	: a < 0 && b < 0 ? 1
147	: -1;
148	}
149
150	template <typename T>
151	static inline
152	int sign_of_compare(T const& a, T const& b, T const& c, T const& d)
153	{
154	// Both ab and cd are positive
155	// We have to judge if ab > cd
156
157	using side::detail::multiplicable_integral;
158	multiplicable_integral<T> ab = multiplicable_integral<T>(a)
159	* multiplicable_integral<T>(b);
160	multiplicable_integral<T> cd = multiplicable_integral<T>(c)
161	* multiplicable_integral<T>(d);
162
163	int result = ab > cd ? 1
164	: ab < cd ? -1
165	: 0
166	;
167
168	#ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
169	using namespace boost::multiprecision;
170	cpp_int const lab = cpp_int(a) * cpp_int(b);
171	cpp_int const lcd = cpp_int(c) * cpp_int(d);
172
173	ab.check_value(lab);
174	cd.check_value(lcd);
175
176	int result2 = lab > lcd ? 1
177	: lab < lcd ? -1
178	: 0
179	;
180	BOOST_GEOMETRY_ASSERT(result == result2);
181	#endif
182
183	return result;
184	}
185
186	template <typename T>
187	static inline
188	int sign_of_addition_of_two_products(T const& a, T const& b, T const& c, T const& d)
189	{
190	// sign of ab+cd, 1 if positive, -1 if negative, else 0
191	int const ab = sign_of_product(a, b);
192	int const cd = sign_of_product(c, d);
193	if (ab == 0)
194	{
195	return cd;
196	}
197	if (cd == 0)
198	{
199	return ab;
200	}
201
202	if (ab == cd)
203	{
204	// Both positive or both negative
205	return ab;
206	}
207
208	// One is positive, one is negative, both are non zero
209	// If ab is positive, we have to judge if ab > -cd (then 1 because sum is positive)
210	// If ab is negative, we have to judge if cd > -ab (idem)
211	return ab == 1
212	? sign_of_compare(a, b, -c, d)
213	: sign_of_compare(c, d, -a, b);
214	}
215
216
217	public :
218
219	// Calculates the side of the intersection-point (if any) of
220	// of segment a//b w.r.t. segment c
221	// This is calculated without (re)calculating the IP itself again and fully
222	// based on integer mathematics
223	template <typename T, typename Segment, typename Point>
224	static inline T side_value(Segment const& a, Segment const& b,
225	Segment const& c, Point const& fallback_point)
226	{
227	// The first point of the three segments is reused several times
228	T const ax = get<0, 0>(a);
229	T const ay = get<0, 1>(a);
230	T const bx = get<0, 0>(b);
231	T const by = get<0, 1>(b);
232	T const cx = get<0, 0>(c);
233	T const cy = get<0, 1>(c);
234
235	T const dx_a = get<1, 0>(a) - ax;
236	T const dy_a = get<1, 1>(a) - ay;
237
238	T const dx_b = get<1, 0>(b) - bx;
239	T const dy_b = get<1, 1>(b) - by;
240
241	T const dx_c = get<1, 0>(c) - cx;
242	T const dy_c = get<1, 1>(c) - cy;
243
244	// Cramer's rule: d (see cart_intersect.hpp)
245	T const d = geometry::detail::determinant<T>
246	(
247	dx_a, dy_a,
248	dx_b, dy_b
249	);
250
251	T const zero = T();
252	if (d == zero)
253	{
254	// There is no IP of a//b, they are collinear or parallel
255	// Assuming they intersect (this method should be called for
256	// segments known to intersect), they are collinear and overlap.
257	// They have one or two intersection points - we don't know and
258	// have to rely on the fallback intersection point
259
260	Point c1, c2;
261	geometry::detail::assign_point_from_index<0>(c, c1);
262	geometry::detail::assign_point_from_index<1>(c, c2);
263	return side_by_triangle<>::apply(c1, c2, fallback_point);
264	}
265
266	// Cramer's rule: da (see cart_intersect.hpp)
267	T const da = geometry::detail::determinant<T>
268	(
269	dx_b, dy_b,
270	ax - bx, ay - by
271	);
272
273	// IP is at (ax + (da/d) * dx_a, ay + (da/d) * dy_a)
274	// Side of IP is w.r.t. c is: determinant(dx_c, dy_c, ipx-cx, ipy-cy)
275	// We replace ipx by expression above and multiply each term by d
276
277	#ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
278	T const result1 = geometry::detail::determinant<T>
279	(
280	dx_c * d, dy_c * d,
281	d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da
282	);
283
284	// Note: result / (d * d)
285	// is identical to the side_value of side_by_triangle
286	// Therefore, the sign is always the same as that result, and the
287	// resulting side (left,right,collinear) is the same
288
289	// The first row we divide again by d because of determinant multiply rule
290	T const result2 = d * geometry::detail::determinant<T>
291	(
292	dx_c, dy_c,
293	d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da
294	);
295	// Write out:
296	T const result3 = d * (dx_c * (d * (ay - cy) + dy_a * da)
297	- dy_c * (d * (ax - cx) + dx_a * da));
298	// Write out in braces:
299	T const result4 = d * (dx_c * d * (ay - cy) + dx_c * dy_a * da
300	- dy_c * d * (ax - cx) - dy_c * dx_a * da);
301	// Write in terms of d * XX + da * YY
302	T const result5 = d * (d * (dx_c * (ay - cy) - dy_c * (ax - cx))
303	+ da * (dx_c * dy_a - dy_c * dx_a));
304
305	boost::ignore_unused(result1, result2, result3, result4, result5);
306	//return result;
307	#endif
308
309	// We consider the results separately
310	// (in the end we only have to return the side-value 1,0 or -1)
311
312	// To avoid multiplications we judge the product (easy, avoids *d)
313	// and the sign of pq+rs (more elaborate)
314	T const result = sign_of_product
315	(
316	d,
317	sign_of_addition_of_two_products
318	(
319	d, dx_c * (ay - cy) - dy_c * (ax - cx),
320	da, dx_c * dy_a - dy_c * dx_a
321	)
322	);
323	return result;
324
325
326	}
327
328	template <typename Segment, typename Point>
329	static inline int apply(Segment const& a, Segment const& b,
330	Segment const& c,
331	Point const& fallback_point)
332	{
333	typedef typename geometry::coordinate_type<Segment>::type coordinate_type;
334	coordinate_type const s = side_value<coordinate_type>(a, b, c, fallback_point);
335	coordinate_type const zero = coordinate_type();
336	return math::equals(s, zero) ? 0
337	: s > zero ? 1
338	: -1;
339	}
340
341	};
342
343
344	}} // namespace strategy::side
345
346	}} // namespace boost::geometry
347
348
349	#endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP