• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 
stopwatchstopwatch64    stopwatch() : m_start(clock_type::now()) { }
65 
stopwatchstopwatch66    stopwatch(const stopwatch& other) : m_start(other.m_start) { }
67 
operator =stopwatch68    stopwatch& operator=(const stopwatch& other)
69    {
70       m_start = other.m_start;
71       return *this;
72    }
73 
~stopwatchstopwatch74    ~stopwatch() { }
75 
elapsedstopwatch76    float elapsed() const
77    {
78       return float(std::chrono::nanoseconds((clock_type::now() - m_start)).count()) / 1e9f;
79    }
80 
resetstopwatch81    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;
scoped_precisionscoped_precision93    scoped_precision(unsigned new_p) : p(mpfr_float::default_precision())
94    {
95       mpfr_float::default_precision(new_p);
96    }
~scoped_precisionscoped_precision97    ~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 */
beta_distribution_median_method_1(mpfr_float const & a_,mpfr_float const & b_,unsigned digits10)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 crude: 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 its derivative in one go:
141 
142 */
beta_distribution_median_method_2(mpfr_float const & a_,mpfr_float const & b_,unsigned digits10)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 methods 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 
beta_distribution_median_method_3(mpfr_float const & a_,mpfr_float const & b_,unsigned digits10)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 
main()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 //] //[/mpfr_variable]
252 
253