[ceph.git] / ceph / src / boost / 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, 2019.
// Modifications copyright (c) 2015-2019, 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;
    }

#ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
    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);
    }
#endif // BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
};

} // 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
92f5a8d4 TL	5	// This file was modified by Oracle on 2015, 2019.
92f5a8d4 TL	6	// Modifications copyright (c) 2015-2019, Oracle and/or its affiliates.
7c673cae FG	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
92f5a8d4	121	#ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
7c673cae FG	122	template <typename CmpVal>
	123	void check_value(CmpVal const& cmp_val) const
	124	{
	125	unsigned_type b = base; // a workaround for MinGW - undefined reference base
	126	CmpVal val = CmpVal(m_sign) * (CmpVal(m_ms) * CmpVal(b) + CmpVal(m_ls));
	127	BOOST_GEOMETRY_ASSERT(cmp_val == val);
	128	}
92f5a8d4	129	#endif // BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
7c673cae FG	130	};
	131
	132	} // namespace detail
	133
	134	// Calculates the side of the intersection-point (if any) of
	135	// of segment a//b w.r.t. segment c
	136	// This is calculated without (re)calculating the IP itself again and fully
	137	// based on integer mathematics; there are no divisions
	138	// It can be used for either integer (rescaled) points, and also for FP
	139	class side_of_intersection
	140	{
	141	private :
	142	template <typename T, typename U>
	143	static inline
	144	int sign_of_product(T const& a, U const& b)
	145	{
	146	return a == 0 \|\| b == 0 ? 0
	147	: a > 0 && b > 0 ? 1
	148	: a < 0 && b < 0 ? 1
	149	: -1;
	150	}
	151
	152	template <typename T>
	153	static inline
	154	int sign_of_compare(T const& a, T const& b, T const& c, T const& d)
	155	{
	156	// Both ab and cd are positive
	157	// We have to judge if ab > cd
	158
	159	using side::detail::multiplicable_integral;
	160	multiplicable_integral<T> ab = multiplicable_integral<T>(a)
	161	* multiplicable_integral<T>(b);
	162	multiplicable_integral<T> cd = multiplicable_integral<T>(c)
	163	* multiplicable_integral<T>(d);
	164
	165	int result = ab > cd ? 1
	166	: ab < cd ? -1
	167	: 0
	168	;
	169
	170	#ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
	171	using namespace boost::multiprecision;
	172	cpp_int const lab = cpp_int(a) * cpp_int(b);
	173	cpp_int const lcd = cpp_int(c) * cpp_int(d);
	174
	175	ab.check_value(lab);
	176	cd.check_value(lcd);
	177
	178	int result2 = lab > lcd ? 1
	179	: lab < lcd ? -1
	180	: 0
	181	;
	182	BOOST_GEOMETRY_ASSERT(result == result2);
	183	#endif
	184
	185	return result;
	186	}
	187
	188	template <typename T>
	189	static inline
	190	int sign_of_addition_of_two_products(T const& a, T const& b, T const& c, T const& d)
	191	{
	192	// sign of ab+cd, 1 if positive, -1 if negative, else 0
	193	int const ab = sign_of_product(a, b);
194	int const cd = sign_of_product(c, d);
195	if (ab == 0)
196	{
197	return cd;
198	}
199	if (cd == 0)
200	{
201	return ab;
202	}
203
204	if (ab == cd)
205	{
206	// Both positive or both negative
207	return ab;
208	}
209
210	// One is positive, one is negative, both are non zero
211	// If ab is positive, we have to judge if ab > -cd (then 1 because sum is positive)
212	// If ab is negative, we have to judge if cd > -ab (idem)
213	return ab == 1
214	? sign_of_compare(a, b, -c, d)
215	: sign_of_compare(c, d, -a, b);
216	}
217
218
219	public :
220
221	// Calculates the side of the intersection-point (if any) of
222	// of segment a//b w.r.t. segment c
223	// This is calculated without (re)calculating the IP itself again and fully
224	// based on integer mathematics
225	template <typename T, typename Segment, typename Point>
226	static inline T side_value(Segment const& a, Segment const& b,
227	Segment const& c, Point const& fallback_point)
228	{
229	// The first point of the three segments is reused several times
230	T const ax = get<0, 0>(a);
231	T const ay = get<0, 1>(a);
232	T const bx = get<0, 0>(b);
233	T const by = get<0, 1>(b);
234	T const cx = get<0, 0>(c);
235	T const cy = get<0, 1>(c);
236
237	T const dx_a = get<1, 0>(a) - ax;
238	T const dy_a = get<1, 1>(a) - ay;
239
240	T const dx_b = get<1, 0>(b) - bx;
241	T const dy_b = get<1, 1>(b) - by;
242
243	T const dx_c = get<1, 0>(c) - cx;
244	T const dy_c = get<1, 1>(c) - cy;
245
246	// Cramer's rule: d (see cart_intersect.hpp)
247	T const d = geometry::detail::determinant<T>
248	(
249	dx_a, dy_a,
250	dx_b, dy_b
251	);
252
253	T const zero = T();
254	if (d == zero)
255	{
256	// There is no IP of a//b, they are collinear or parallel
257	// Assuming they intersect (this method should be called for
258	// segments known to intersect), they are collinear and overlap.
259	// They have one or two intersection points - we don't know and
260	// have to rely on the fallback intersection point
261
262	Point c1, c2;
263	geometry::detail::assign_point_from_index<0>(c, c1);
264	geometry::detail::assign_point_from_index<1>(c, c2);
265	return side_by_triangle<>::apply(c1, c2, fallback_point);
266	}
267
268	// Cramer's rule: da (see cart_intersect.hpp)
269	T const da = geometry::detail::determinant<T>
270	(
271	dx_b, dy_b,
272	ax - bx, ay - by
273	);
274
275	// IP is at (ax + (da/d) * dx_a, ay + (da/d) * dy_a)
276	// Side of IP is w.r.t. c is: determinant(dx_c, dy_c, ipx-cx, ipy-cy)
277	// We replace ipx by expression above and multiply each term by d
278
279	#ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG
280	T const result1 = geometry::detail::determinant<T>
281	(
282	dx_c * d, dy_c * d,
283	d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da
284	);
285
286	// Note: result / (d * d)
287	// is identical to the side_value of side_by_triangle
288	// Therefore, the sign is always the same as that result, and the
289	// resulting side (left,right,collinear) is the same
290
291	// The first row we divide again by d because of determinant multiply rule
292	T const result2 = d * geometry::detail::determinant<T>
293	(
294	dx_c, dy_c,
295	d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da
296	);
297	// Write out:
298	T const result3 = d * (dx_c * (d * (ay - cy) + dy_a * da)
299	- dy_c * (d * (ax - cx) + dx_a * da));
300	// Write out in braces:
301	T const result4 = d * (dx_c * d * (ay - cy) + dx_c * dy_a * da
302	- dy_c * d * (ax - cx) - dy_c * dx_a * da);
303	// Write in terms of d * XX + da * YY
304	T const result5 = d * (d * (dx_c * (ay - cy) - dy_c * (ax - cx))
305	+ da * (dx_c * dy_a - dy_c * dx_a));
306
307	boost::ignore_unused(result1, result2, result3, result4, result5);
308	//return result;
309	#endif
310
311	// We consider the results separately
312	// (in the end we only have to return the side-value 1,0 or -1)
313
314	// To avoid multiplications we judge the product (easy, avoids *d)
315	// and the sign of pq+rs (more elaborate)
316	T const result = sign_of_product
317	(
318	d,
319	sign_of_addition_of_two_products
320	(
321	d, dx_c * (ay - cy) - dy_c * (ax - cx),
322	da, dx_c * dy_a - dy_c * dx_a
323	)
324	);
325	return result;
326
327
328	}
329
330	template <typename Segment, typename Point>
331	static inline int apply(Segment const& a, Segment const& b,
332	Segment const& c,
333	Point const& fallback_point)
334	{
335	typedef typename geometry::coordinate_type<Segment>::type coordinate_type;
336	coordinate_type const s = side_value<coordinate_type>(a, b, c, fallback_point);
337	coordinate_type const zero = coordinate_type();
338	return math::equals(s, zero) ? 0
339	: s > zero ? 1
340	: -1;
341	}
342
343	};
344
345
346	}} // namespace strategy::side
347
348	}} // namespace boost::geometry
349
350
351	#endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP