]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
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 | |
f67539c2 | 138 | crude: we simply assume that we have approximately 14 digits correct from the double-precision approximation and |
92f5a8d4 | 139 | that the precision doubles with each step. We also cheat, and use an internal Boost.Math function that calculates |
f67539c2 | 140 | ['I[sub x](a, b)] and its derivative in one go: |
92f5a8d4 TL |
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 | /*` | |
f67539c2 | 159 | Before we refine the method further, it might be wise to take stock and see how methods 1 and 2 compare. |
92f5a8d4 TL |
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 | } | |
f67539c2 | 251 | //] //[/mpfr_variable] |
92f5a8d4 | 252 |