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