libs/multiprecision/performance/delaunay_test.cpp

   1 ///////////////////////////////////////////////////////////////////////////////
   2 //  Copyright 2012 John Maddock.
   3 //  Copyright 2012 Phil Endecott
   4 //  Distributed under the Boost
   5 //  Software License, Version 1.0. (See accompanying file
   6 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   7
   8 #include <boost/multiprecision/cpp_int.hpp>
   9 #include "arithmetic_backend.hpp"
  10 #include <boost/chrono.hpp>
  11 #include <boost/random/mersenne_twister.hpp>
  12 #include <boost/random/uniform_int_distribution.hpp>
  13
  14 #include <fstream>
  15 #include <iomanip>
  16
  17 template <class Clock>
  18 struct stopwatch
  19 {
  20    typedef typename Clock::duration duration;
  21    stopwatch()
  22    {
  23       m_start = Clock::now();
  24    }
  25    duration elapsed()
  26    {
  27       return Clock::now() - m_start;
  28    }
  29    void reset()
  30    {
  31       m_start = Clock::now();
  32    }
  33
  34  private:
  35    typename Clock::time_point m_start;
  36 };
  37
  38 // Custom 128-bit maths used for exact calculation of the Delaunay test.
  39 // Only the few operators actually needed here are implemented.
  40
  41 struct int128_t
  42 {
  43    int64_t  high;
  44    uint64_t low;
  45
  46    int128_t() {}
  47    int128_t(int32_t i) : high(i >> 31), low(static_cast<int64_t>(i)) {}
  48    int128_t(uint32_t i) : high(0), low(i) {}
  49    int128_t(int64_t i) : high(i >> 63), low(i) {}
  50    int128_t(uint64_t i) : high(0), low(i) {}
  51 };
  52
  53 inline int128_t operator<<(int128_t val, int amt)
  54 {
  55    int128_t r;
  56    r.low  = val.low << amt;
  57    r.high = val.low >> (64 - amt);
  58    r.high |= val.high << amt;
  59    return r;
  60 }
  61
  62 inline int128_t& operator+=(int128_t& l, int128_t r)
  63 {
  64    l.low += r.low;
  65    bool carry = l.low < r.low;
  66    l.high += r.high;
  67    if (carry)
  68       ++l.high;
  69    return l;
  70 }
  71
  72 inline int128_t operator-(int128_t val)
  73 {
  74    val.low  = ~val.low;
  75    val.high = ~val.high;
  76    val.low += 1;
  77    if (val.low == 0)
  78       val.high += 1;
  79    return val;
  80 }
  81
  82 inline int128_t operator+(int128_t l, int128_t r)
  83 {
  84    l += r;
  85    return l;
  86 }
  87
  88 inline bool operator<(int128_t l, int128_t r)
  89 {
  90    if (l.high != r.high)
  91       return l.high < r.high;
  92    return l.low < r.low;
  93 }
  94
  95 inline int128_t mult_64x64_to_128(int64_t a, int64_t b)
  96 {
  97    // Make life simple by dealing only with positive numbers:
  98    bool neg = false;
  99    if (a < 0)
 100    {
 101       neg = !neg;
 102       a   = -a;
 103    }
 104    if (b < 0)
 105    {
 106       neg = !neg;
 107       b   = -b;
 108    }
 109
 110    // Divide input into 32-bit halves:
 111    uint32_t ah = a >> 32;
 112    uint32_t al = a & 0xffffffff;
 113    uint32_t bh = b >> 32;
 114    uint32_t bl = b & 0xffffffff;
 115
 116    // Long multiplication, with 64-bit temporaries:
 117
 118    //            ah al
 119    //          * bh bl
 120    // ----------------
 121    //            al*bl   (t1)
 122    // +       ah*bl      (t2)
 123    // +       al*bh      (t3)
 124    // +    ah*bh         (t4)
 125    // ----------------
 126
 127    uint64_t t1 = static_cast<uint64_t>(al) * bl;
 128    uint64_t t2 = static_cast<uint64_t>(ah) * bl;
 129    uint64_t t3 = static_cast<uint64_t>(al) * bh;
 130    uint64_t t4 = static_cast<uint64_t>(ah) * bh;
 131
 132    int128_t r(t1);
 133    r.high = t4;
 134    r += int128_t(t2) << 32;
 135    r += int128_t(t3) << 32;
 136
 137    if (neg)
 138       r = -r;
 139
 140    return r;
 141 }
 142
 143 template <class R, class T>
 144 BOOST_FORCEINLINE void mul_2n(R& r, const T& a, const T& b)
 145 {
 146    r = a;
 147    r *= b;
 148 }
 149
 150 template <class B, boost::multiprecision::expression_template_option ET, class T>
 151 BOOST_FORCEINLINE void mul_2n(boost::multiprecision::number<B, ET>& r, const T& a, const T& b)
 152 {
 153    multiply(r, a, b);
 154 }
 155
 156 BOOST_FORCEINLINE void mul_2n(int128_t& r, const boost::int64_t& a, const boost::int64_t& b)
 157 {
 158    r = mult_64x64_to_128(a, b);
 159 }
 160
 161 template <class Traits>
 162 inline bool delaunay_test(int32_t ax, int32_t ay, int32_t bx, int32_t by,
 163                           int32_t cx, int32_t cy, int32_t dx, int32_t dy)
 164 {
 165    // Test whether the quadrilateral ABCD's diagonal AC should be flipped to BD.
 166    // This is the Cline & Renka method.
 167    // Flip if the sum of the angles ABC and CDA is greater than 180 degrees.
 168    // Equivalently, flip if sin(ABC + CDA) < 0.
 169    // Trig identity: cos(ABC) * sin(CDA) + sin(ABC) * cos(CDA) < 0
 170    // We can use scalar and vector products to find sin and cos, and simplify
 171    // to the following code.
 172    // Numerical robustness is important.  This code addresses it by performing
 173    // exact calculations with large integer types.
 174    //
 175    // NOTE: This routine is limited to inputs with up to 30 BIT PRECISION, which
 176    // is to say all inputs must be in the range [INT_MIN/2, INT_MAX/2].
 177
 178    typedef typename Traits::i64_t  i64;
 179    typedef typename Traits::i128_t i128;
 180
 181    i64 cos_abc, t;
 182    mul_2n(cos_abc, (ax - bx), (cx - bx)); // subtraction yields 31-bit values, multiplied to give 62-bit values
 183    mul_2n(t, (ay - by), (cy - by));
 184    cos_abc += t; // addition yields 63 bit value, leaving one left for the sign
 185
 186    i64 cos_cda;
 187    mul_2n(cos_cda, (cx - dx), (ax - dx));
 188    mul_2n(t, (cy - dy), (ay - dy));
 189    cos_cda += t;
 190
 191    if (cos_abc >= 0 && cos_cda >= 0)
 192       return false;
 193    if (cos_abc < 0 && cos_cda < 0)
 194       return true;
 195
 196    i64 sin_abc;
 197    mul_2n(sin_abc, (ax - bx), (cy - by));
 198    mul_2n(t, (cx - bx), (ay - by));
 199    sin_abc -= t;
 200
 201    i64 sin_cda;
 202    mul_2n(sin_cda, (cx - dx), (ay - dy));
 203    mul_2n(t, (ax - dx), (cy - dy));
 204    sin_cda -= t;
 205
 206    i128 sin_sum, t128;
 207    mul_2n(sin_sum, sin_abc, cos_cda); // 63-bit inputs multiplied to 126-bit output
 208    mul_2n(t128, cos_abc, sin_cda);
 209    sin_sum += t128; // Addition yields 127 bit result, leaving one bit for the sign
 210
 211    return sin_sum < 0;
 212 }
 213
 214 struct dt_dat
 215 {
 216    int32_t ax, ay, bx, by, cx, cy, dx, dy;
 217 };
 218
 219 typedef std::vector<dt_dat> data_t;
 220 data_t                      data;
 221
 222 template <class Traits>
 223 void do_calc(const char* name)
 224 {
 225    std::cout << "Running calculations for: " << name << std::endl;
 226
 227    stopwatch<boost::chrono::high_resolution_clock> w;
 228
 229    boost::uint64_t flips = 0;
 230    boost::uint64_t calcs = 0;
 231
 232    for (int j = 0; j < 1000; ++j)
 233    {
 234       for (data_t::const_iterator i = data.begin(); i != data.end(); ++i)
 235       {
 236          const dt_dat& d    = *i;
 237          bool          flip = delaunay_test<Traits>(d.ax, d.ay, d.bx, d.by, d.cx, d.cy, d.dx, d.dy);
 238          if (flip)
 239             ++flips;
 240          ++calcs;
 241       }
 242    }
 243    double t = boost::chrono::duration_cast<boost::chrono::duration<double> >(w.elapsed()).count();
 244
 245    std::cout << "Number of calculations = " << calcs << std::endl;
 246    std::cout << "Number of flips = " << flips << std::endl;
 247    std::cout << "Total execution time = " << t << std::endl;
 248    std::cout << "Time per calculation = " << t / calcs << std::endl
 249              << std::endl;
 250 }
 251
 252 template <class I64, class I128>
 253 struct test_traits
 254 {
 255    typedef I64  i64_t;
 256    typedef I128 i128_t;
 257 };
 258
 259 dt_dat generate_quadrilateral()
 260 {
 261    static boost::random::mt19937                    gen;
 262    static boost::random::uniform_int_distribution<> dist(INT_MIN / 2, INT_MAX / 2);
 263
 264    dt_dat result;
 265
 266    result.ax = dist(gen);
 267    result.ay = dist(gen);
 268    result.bx = boost::random::uniform_int_distribution<>(result.ax, INT_MAX / 2)(gen); // bx is to the right of ax.
 269    result.by = dist(gen);
 270    result.cx = dist(gen);
 271    result.cy = boost::random::uniform_int_distribution<>(result.cx > result.bx ? result.by : result.ay, INT_MAX / 2)(gen); // cy is below at least one of ay and by.
 272    result.dx = boost::random::uniform_int_distribution<>(result.cx, INT_MAX / 2)(gen);                                     // dx is to the right of cx.
 273    result.dy = boost::random::uniform_int_distribution<>(result.cx > result.bx ? result.by : result.ay, INT_MAX / 2)(gen); // cy is below at least one of ay and by.
 274
 275    return result;
 276 }
 277
 278 static void load_data()
 279 {
 280    for (unsigned i = 0; i < 100000; ++i)
 281       data.push_back(generate_quadrilateral());
 282 }
 283
 284 int main()
 285 {
 286    using namespace boost::multiprecision;
 287    std::cout << "loading data...\n";
 288    load_data();
 289
 290    std::cout << "calculating...\n";
 291
 292    do_calc<test_traits<boost::int64_t, boost::int64_t> >("int64_t, int64_t");
 293    do_calc<test_traits<number<arithmetic_backend<boost::int64_t>, et_off>, number<arithmetic_backend<boost::int64_t>, et_off> > >("arithmetic_backend<int64_t>, arithmetic_backend<int64_t>");
 294    do_calc<test_traits<boost::int64_t, number<arithmetic_backend<boost::int64_t>, et_off> > >("int64_t, arithmetic_backend<int64_t>");
 295    do_calc<test_traits<number<cpp_int_backend<64, 64, boost::multiprecision::signed_magnitude, boost::multiprecision::unchecked, void>, et_off>, number<cpp_int_backend<64, 64, boost::multiprecision::signed_magnitude, boost::multiprecision::unchecked, void>, et_off> > >("multiprecision::int64_t, multiprecision::int64_t");
 296
 297    do_calc<test_traits<boost::int64_t, ::int128_t> >("int64_t, int128_t");
 298    do_calc<test_traits<boost::int64_t, boost::multiprecision::int128_t> >("int64_t, boost::multiprecision::int128_t");
 299    do_calc<test_traits<boost::int64_t, number<cpp_int_backend<128, 128, boost::multiprecision::signed_magnitude, boost::multiprecision::unchecked, void>, et_on> > >("int64_t, int128_t (ET)");
 300    do_calc<test_traits<number<cpp_int_backend<64, 64, boost::multiprecision::signed_magnitude, boost::multiprecision::unchecked, void>, et_off>, boost::multiprecision::int128_t> >("multiprecision::int64_t, multiprecision::int128_t");
 301
 302    do_calc<test_traits<boost::int64_t, cpp_int> >("int64_t, cpp_int");
 303    do_calc<test_traits<boost::int64_t, number<cpp_int_backend<>, et_off> > >("int64_t, cpp_int (no ET's)");
 304    do_calc<test_traits<boost::int64_t, number<cpp_int_backend<128> > > >("int64_t, cpp_int(128-bit cache)");
 305    do_calc<test_traits<boost::int64_t, number<cpp_int_backend<128>, et_off> > >("int64_t, cpp_int (128-bit Cache no ET's)");
 306
 307    return 0;
 308 }