libs/multiprecision/example/mpfr_precision.cpp

   1 ///////////////////////////////////////////////////////////////
   2 //  Copyright 2018 John Maddock. Distributed under the Boost
   3 //  Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
   5
   6 //[mpfr_variable
   7
   8 /*`
   9 This example illustrates the use of variable-precision arithmetic with
  10 the `mpfr_float` number type.  We'll calculate the median of the
  11 beta distribution to an absurdly high precision and compare the
  12 accuracy and times taken for various methods.  That is, we want
  13 to calculate the value of `x` for which ['I[sub x](a, b) = 0.5].
  14
  15 Ultimately we'll use Newtons method and set the precision of
  16 mpfr_float to have just enough digits at each iteration.
  17
  18 The full source of the this program is in [@../../example/mpfr_precision.cpp]
  19
  20 We'll skip over the #includes and using declations, and go straight to
  21 some support code, first off a simple stopwatch for performance measurement:
  22
  23 */
  24
  25 //=template <class clock_type>
  26 //=struct stopwatch { /*details \*/ };
  27
  28 /*`
  29 We'll use `stopwatch<std::chono::high_resolution_clock>` as our performance measuring device.
  30
  31 We also have a small utility class for controlling the current precision of mpfr_float:
  32
  33    struct scoped_precision
  34    {
  35       unsigned p;
  36       scoped_precision(unsigned new_p) : p(mpfr_float::default_precision())
  37       {
  38          mpfr_float::default_precision(new_p);
  39       }
  40       ~scoped_precision()
  41       {
  42          mpfr_float::default_precision(p);
  43       }
  44    };
  45
  46 */
  47 //<-
  48 #include <boost/multiprecision/mpfr.hpp>
  49 #include <boost/math/special_functions/beta.hpp>
  50 #include <boost/math/special_functions/relative_difference.hpp>
  51 #include <iostream>
  52 #include <chrono>
  53
  54 using boost::multiprecision::mpfr_float;
  55 using boost::math::ibeta_inv;
  56 using namespace boost::math::policies;
  57
  58 template <class clock_type>
  59 struct stopwatch
  60 {
  61 public:
  62    typedef typename clock_type::duration duration_type;
  63
  64    stopwatch() : m_start(clock_type::now()) { }
  65
  66    stopwatch(const stopwatch& other) : m_start(other.m_start) { }
  67
  68    stopwatch& operator=(const stopwatch& other)
  69    {
  70       m_start = other.m_start;
  71       return *this;
  72    }
  73
  74    ~stopwatch() { }
  75
  76    float elapsed() const
  77    {
  78       return float(std::chrono::nanoseconds((clock_type::now() - m_start)).count()) / 1e9f;
  79    }
  80
  81    void reset()
  82    {
  83       m_start = clock_type::now();
  84    }
  85
  86 private:
  87    typename clock_type::time_point m_start;
  88 };
  89
  90 struct scoped_precision
  91 {
  92    unsigned p;
  93    scoped_precision(unsigned new_p) : p(mpfr_float::default_precision())
  94    {
  95       mpfr_float::default_precision(new_p);
  96    }
  97    ~scoped_precision()
  98    {
  99       mpfr_float::default_precision(p);
 100    }
 101 };
 102 //->
 103
 104 /*`
 105 We'll begin with a reference method that simply calls the Boost.Math function `ibeta_inv` and uses the
 106 full working precision of the arguments throughout.  Our reference function takes 3 arguments:
 107
 108 * The 2 parameters `a` and `b` of the beta distribution, and
 109 * The number of decimal digits precision to achieve in the result.
 110
 111 We begin by setting the default working precision to that requested, and then, since we don't know where
 112 our arguments `a` and `b` have been or what precision they have, we make a copy of them - note that since
 113 copying also copies the precision as well as the value, we have to set the precision expicitly with a
 114 second argument to the copy.  Then we can simply return the result of `ibeta_inv`:
 115 */
 116 mpfr_float beta_distribution_median_method_1(mpfr_float const& a_, mpfr_float const& b_, unsigned digits10)
 117 {
 118    scoped_precision sp(digits10);
 119    mpfr_float half(0.5), a(a_, digits10), b(b_, digits10);
 120    return ibeta_inv(a, b, half);
 121 }
 122 /*`
 123 You be wondering why we needed to change the precision of our variables `a` and `b` as well as setting the default -
 124 there are in fact two ways in which this can go wrong if we don't do that:
 125
 126 * The variables have too much precision - this will cause all arithmetic operations involving those types to be
 127 promoted to the higher precision wasting precious calculation time.
 128 * The variables have too little precision - this will cause expressions involving only those variables to be
 129 calculated at the lower precision - for example if we calculate `exp(a)` internally, this will be evaluated at
 130 the precision of `a`, and not the current default.
 131
 132 Since our reference method carries out all calculations at the full precision requested, an obvious refinement
 133 would be to calculate a first approximation to `double` precision and then to use Newton steps to refine it further.
 134
 135 Our function begins the same as before: set the new default precision and then make copies of our arguments
 136 at the correct precision.  We then call `ibeta_inv` with all double precision arguments, promote the result
 137 to an `mpfr_float` and perform Newton steps to obtain the result.  Note that our termination condition is somewhat
 138 cude: we simply assume that we have approximately 14 digits correct from the double-precision approximation and
 139 that the precision doubles with each step.  We also cheat, and use an internal Boost.Math function that calculates
 140 ['I[sub x](a, b)] and it's derivative in one go:
 141
 142 */
 143 mpfr_float beta_distribution_median_method_2(mpfr_float const& a_, mpfr_float const& b_, unsigned digits10)
 144 {
 145    scoped_precision sp(digits10);
 146    mpfr_float half(0.5), a(a_, digits10), b(b_, digits10);
 147    mpfr_float guess = ibeta_inv((double)a, (double)b, 0.5);
 148    unsigned current_digits = 14;
 149    mpfr_float f, f1;
 150    while (current_digits < digits10)
 151    {
 152       f = boost::math::detail::ibeta_imp(a, b, guess, boost::math::policies::policy<>(), false, true, &f1) - half;
 153       guess -= f / f1;
 154       current_digits *= 2;
 155    }
 156    return guess;
 157 }
 158 /*`
 159 Before we refine the method further, it might be wise to take stock and see how method's 1 and 2 compare.
 160 We'll ask them both for 1500 digit precision, and compare against the value produced by `ibeta_inv` at 1700 digits.
 161 Here's what the results look like:
 162
 163 [pre
 164 Method 1 time = 0.611647
 165 Relative error: 2.99991e-1501
 166 Method 2 time = 0.646746
 167 Relative error: 7.55843e-1501
 168 ]
 169
 170 Clearly they are both equally accurate, but Method 1 is actually faster and our plan for improved performance
 171 hasn't actually worked.  It turns out that we're not actually comparing like with like, because `ibeta_inv` uses
 172 Halley iteration internally which churns out more digits of precision rather more rapidly than Newton iteration.
 173 So the time we save by refining an initial `double` approximation, then loose it again by taking more iterations
 174 to get to the result.
 175
 176 Time for a more refined approach.  It follows the same form as Method 2, but now we set the working precision
 177 within the Newton iteration loop, to just enough digits to cover the expected precision at each step.  That means
 178 we also create new copies of our arguments at the correct precision within the loop, and likewise change the precision
 179 of the current `guess` each time through:
 180
 181 */
 182
 183 mpfr_float beta_distribution_median_method_3(mpfr_float const& a_, mpfr_float const& b_, unsigned digits10)
 184 {
 185    mpfr_float guess = ibeta_inv((double)a_, (double)b_, 0.5);
 186    unsigned current_digits = 14;
 187    mpfr_float f(0, current_digits), f1(0, current_digits), delta(1);
 188    while (current_digits < digits10)
 189    {
 190       current_digits *= 2;
 191       scoped_precision sp((std::min)(current_digits, digits10));
 192       mpfr_float a(a_, mpfr_float::default_precision()), b(b_, mpfr_float::default_precision());
 193       guess.precision(mpfr_float::default_precision());
 194       f = boost::math::detail::ibeta_imp(a, b, guess, boost::math::policies::policy<>(), false, true, &f1) - 0.5f;
 195       guess -= f / f1;
 196    }
 197    return guess;
 198 }
 199
 200 /*`
 201 The new performance results look much more promising:
 202
 203 [pre
 204 Method 1 time = 0.591244
 205 Relative error: 2.99991e-1501
 206 Method 2 time = 0.622679
 207 Relative error: 7.55843e-1501
 208 Method 3 time = 0.143393
 209 Relative error: 4.03898e-1501
 210 ]
 211
 212 This time we're 4x faster than `ibeta_inv`, and no doubt that could be improved a little more by carefully
 213 optimising the number of iterations and the method (Halley vs Newton) taken.
 214
 215 Finally, here's the driver code for the above methods:
 216
 217 */
 218
 219 int main()
 220 {
 221    try {
 222       mpfr_float a(10), b(20);
 223
 224       mpfr_float true_value = beta_distribution_median_method_1(a, b, 1700);
 225
 226       stopwatch<std::chrono::high_resolution_clock> my_stopwatch;
 227
 228       mpfr_float v1 = beta_distribution_median_method_1(a, b, 1500);
 229       float hp_time = my_stopwatch.elapsed();
 230       std::cout << "Method 1 time = " << hp_time << std::endl;
 231       std::cout << "Relative error: " << boost::math::relative_difference(v1, true_value) << std::endl;
 232
 233       my_stopwatch.reset();
 234       mpfr_float v2 = beta_distribution_median_method_2(a, b, 1500);
 235       hp_time = my_stopwatch.elapsed();
 236       std::cout << "Method 2 time = " << hp_time << std::endl;
 237       std::cout << "Relative error: " << boost::math::relative_difference(v2, true_value) << std::endl;
 238
 239       my_stopwatch.reset();
 240       mpfr_float v3 = beta_distribution_median_method_3(a, b, 1500);
 241       hp_time = my_stopwatch.elapsed();
 242       std::cout << "Method 3 time = " << hp_time << std::endl;
 243       std::cout << "Relative error: " << boost::math::relative_difference(v3, true_value) << std::endl;
 244    }
 245    catch (const std::exception& e)
 246    {
 247       std::cout << "Found exception with message: " << e.what() << std::endl;
 248    }
 249    return 0;
 250 }
 251 //]
 252