]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | /////////////////////////////////////////////////////////////// |
2 | // Copyright 2012 John Maddock. Distributed under the Boost | |
3 | // Software License, Version 1.0. (See accompanying file | |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_ | |
5 | ||
6 | #include <boost/math/constants/constants.hpp> | |
7 | #include <boost/multiprecision/cpp_dec_float.hpp> | |
8 | #include <boost/math/special_functions/gamma.hpp> | |
9 | #include <boost/math/special_functions/bessel.hpp> | |
10 | #include <iostream> | |
11 | #include <iomanip> | |
12 | ||
13 | #if !defined(BOOST_NO_CXX11_HDR_ARRAY) && !defined(BOOST_NO_CXX11_LAMBDAS) | |
14 | ||
15 | #include <array> | |
16 | ||
17 | //[AOS1 | |
18 | ||
19 | /*`Generic numeric programming employs templates to use the same code for different | |
20 | floating-point types and functions. Consider the area of a circle a of radius r, given by | |
21 | ||
22 | [:['a = [pi] * r[super 2]]] | |
23 | ||
24 | The area of a circle can be computed in generic programming using Boost.Math | |
25 | for the constant [pi] as shown below: | |
26 | ||
27 | */ | |
28 | ||
29 | //=#include <boost/math/constants/constants.hpp> | |
30 | ||
31 | template<typename T> | |
32 | inline T area_of_a_circle(T r) | |
33 | { | |
34 | using boost::math::constants::pi; | |
35 | return pi<T>() * r * r; | |
36 | } | |
37 | ||
38 | /*` | |
39 | It is possible to use `area_of_a_circle()` with built-in floating-point types as | |
40 | well as floating-point types from Boost.Multiprecision. In particular, consider a | |
41 | system with 4-byte single-precision float, 8-byte double-precision double and also the | |
42 | `cpp_dec_float_50` data type from Boost.Multiprecision with 50 decimal digits | |
43 | of precision. | |
44 | ||
45 | We can compute and print the approximate area of a circle with radius 123/100 for | |
46 | `float`, `double` and `cpp_dec_float_50` with the program below. | |
47 | ||
48 | */ | |
49 | ||
50 | //] | |
51 | ||
52 | //[AOS3 | |
53 | ||
54 | /*`In the next example we'll look at calling both standard library and Boost.Math functions from within generic code. | |
55 | We'll also show how to cope with template arguments which are expression-templates rather than number types.*/ | |
56 | ||
57 | //] | |
58 | ||
59 | //[JEL | |
60 | ||
61 | /*` | |
62 | In this example we'll show several implementations of the | |
63 | [@http://mathworld.wolfram.com/LambdaFunction.html Jahnke and Emden Lambda function], | |
64 | each implementation a little more sophisticated than the last. | |
65 | ||
66 | The Jahnke-Emden Lambda function is defined by the equation: | |
67 | ||
68 | [:['JahnkeEmden(v, z) = [Gamma](v+1) * J[sub v](z) / (z / 2)[super v]]] | |
69 | ||
70 | If we were to implement this at double precision using Boost.Math's facilities for the Gamma and Bessel | |
71 | function calls it would look like this: | |
72 | ||
73 | */ | |
74 | ||
75 | double JEL1(double v, double z) | |
76 | { | |
77 | return boost::math::tgamma(v + 1) * boost::math::cyl_bessel_j(v, z) / std::pow(z / 2, v); | |
78 | } | |
79 | ||
80 | /*` | |
81 | Calling this function as: | |
82 | ||
83 | std::cout << std::scientific << std::setprecision(std::numeric_limits<double>::digits10); | |
84 | std::cout << JEL1(2.5, 0.5) << std::endl; | |
85 | ||
86 | Yields the output: | |
87 | ||
88 | [pre 9.822663964796047e-001] | |
89 | ||
90 | Now let's implement the function again, but this time using the multiprecision type | |
91 | `cpp_dec_float_50` as the argument type: | |
92 | ||
93 | */ | |
94 | ||
95 | boost::multiprecision::cpp_dec_float_50 | |
96 | JEL2(boost::multiprecision::cpp_dec_float_50 v, boost::multiprecision::cpp_dec_float_50 z) | |
97 | { | |
98 | return boost::math::tgamma(v + 1) * boost::math::cyl_bessel_j(v, z) / boost::multiprecision::pow(z / 2, v); | |
99 | } | |
100 | ||
101 | /*` | |
102 | The implementation is almost the same as before, but with one key difference - we can no longer call | |
103 | `std::pow`, instead we must call the version inside the `boost::multiprecision` namespace. In point of | |
104 | fact, we could have omitted the namespace prefix on the call to `pow` since the right overload would | |
105 | have been found via [@http://en.wikipedia.org/wiki/Argument-dependent_name_lookup | |
106 | argument dependent lookup] in any case. | |
107 | ||
108 | Note also that the first argument to `pow` along with the argument to `tgamma` in the above code | |
109 | are actually expression templates. The `pow` and `tgamma` functions will handle these arguments | |
110 | just fine. | |
111 | ||
112 | Here's an example of how the function may be called: | |
113 | ||
114 | std::cout << std::scientific << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10); | |
115 | std::cout << JEL2(cpp_dec_float_50(2.5), cpp_dec_float_50(0.5)) << std::endl; | |
116 | ||
117 | Which outputs: | |
118 | ||
119 | [pre 9.82266396479604757017335009796882833995903762577173e-01] | |
120 | ||
121 | Now that we've seen some non-template examples, lets repeat the code again, but this time as a template | |
122 | that can be called either with a builtin type (`float`, `double` etc), or with a multiprecision type: | |
123 | ||
124 | */ | |
125 | ||
126 | template <class Float> | |
127 | Float JEL3(Float v, Float z) | |
128 | { | |
129 | using std::pow; | |
130 | return boost::math::tgamma(v + 1) * boost::math::cyl_bessel_j(v, z) / pow(z / 2, v); | |
131 | } | |
132 | ||
133 | /*` | |
134 | ||
135 | Once again the code is almost the same as before, but the call to `pow` has changed yet again. | |
136 | We need the call to resolve to either `std::pow` (when the argument is a builtin type), or | |
137 | to `boost::multiprecision::pow` (when the argument is a multiprecision type). We do that by | |
138 | making the call unqualified so that versions of `pow` defined in the same namespace as type | |
139 | `Float` are found via argument dependent lookup, while the `using std::pow` directive makes | |
140 | the standard library versions visible for builtin floating point types. | |
141 | ||
142 | Let's call the function with both `double` and multiprecision arguments: | |
143 | ||
144 | std::cout << std::scientific << std::setprecision(std::numeric_limits<double>::digits10); | |
145 | std::cout << JEL3(2.5, 0.5) << std::endl; | |
146 | std::cout << std::scientific << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10); | |
147 | std::cout << JEL3(cpp_dec_float_50(2.5), cpp_dec_float_50(0.5)) << std::endl; | |
148 | ||
149 | Which outputs: | |
150 | ||
151 | [pre | |
152 | 9.822663964796047e-001 | |
153 | 9.82266396479604757017335009796882833995903762577173e-01 | |
154 | ] | |
155 | ||
156 | Unfortunately there is a problem with this version: if we were to call it like this: | |
157 | ||
158 | boost::multiprecision::cpp_dec_float_50 v(2), z(0.5); | |
159 | JEL3(v + 0.5, z); | |
160 | ||
161 | Then we would get a long and inscrutable error message from the compiler: the problem here is that the first | |
162 | argument to `JEL3` is not a number type, but an expression template. We could obviously add a typecast to | |
163 | fix the issue: | |
164 | ||
165 | JEL(cpp_dec_float_50(v + 0.5), z); | |
166 | ||
167 | However, if we want the function JEL to be truly reusable, then a better solution might be preferred. | |
168 | To achieve this we can borrow some code from Boost.Math which calculates the return type of mixed-argument | |
169 | functions, here's how the new code looks now: | |
170 | ||
171 | */ | |
172 | ||
173 | template <class Float1, class Float2> | |
174 | typename boost::math::tools::promote_args<Float1, Float2>::type | |
175 | JEL4(Float1 v, Float2 z) | |
176 | { | |
177 | using std::pow; | |
178 | return boost::math::tgamma(v + 1) * boost::math::cyl_bessel_j(v, z) / pow(z / 2, v); | |
179 | } | |
180 | ||
181 | /*` | |
182 | ||
183 | As you can see the two arguments to the function are now separate template types, and | |
184 | the return type is computed using the `promote_args` metafunction from Boost.Math. | |
185 | ||
186 | Now we can call: | |
187 | ||
188 | std::cout << std::scientific << std::setprecision(std::numeric_limits<cpp_dec_float_100>::digits10); | |
189 | std::cout << JEL4(cpp_dec_float_100(2) + 0.5, cpp_dec_float_100(0.5)) << std::endl; | |
190 | ||
191 | And get 100 digits of output: | |
192 | ||
193 | [pre 9.8226639647960475701733500979688283399590376257717309069410413822165082248153638454147004236848917775e-01] | |
194 | ||
195 | As a bonus, we can now call the function not just with expression templates, but with other mixed types as well: | |
196 | for example `float` and `double` or `int` and `double`, and the correct return type will be computed in each case. | |
197 | ||
198 | Note that while in this case we didn't have to change the body of the function, in the general case | |
199 | any function like this which creates local variables internally would have to use `promote_args` | |
200 | to work out what type those variables should be, for example: | |
201 | ||
202 | template <class Float1, class Float2> | |
203 | typename boost::math::tools::promote_args<Float1, Float2>::type | |
204 | JEL5(Float1 v, Float2 z) | |
205 | { | |
206 | using std::pow; | |
207 | typedef typename boost::math::tools::promote_args<Float1, Float2>::type variable_type; | |
208 | variable_type t = pow(z / 2, v); | |
209 | return boost::math::tgamma(v + 1) * boost::math::cyl_bessel_j(v, z) / t; | |
210 | } | |
211 | ||
212 | */ | |
213 | ||
214 | //] | |
215 | ||
216 | //[ND1 | |
217 | ||
218 | /*` | |
219 | In this example we'll add even more power to generic numeric programming using not only different | |
220 | floating-point types but also function objects as template parameters. Consider | |
221 | some well-known central difference rules for numerically computing the first derivative | |
222 | of a function ['f[prime](x)] with ['x [isin] [real]]: | |
223 | ||
224 | [equation floating_point_eg1] | |
225 | ||
226 | Where the difference terms ['m[sub n]] are given by: | |
227 | ||
228 | [equation floating_point_eg2] | |
229 | ||
230 | and ['dx] is the step-size of the derivative. | |
231 | ||
232 | The third formula in Equation 1 is a three-point central difference rule. It calculates | |
233 | the first derivative of ['f[prime](x)] to ['O(dx[super 6])], where ['dx] is the given step-size. | |
234 | For example, if | |
235 | the step-size is 0.01 this derivative calculation has about 6 decimal digits of precision - | |
236 | just about right for the 7 decimal digits of single-precision float. | |
237 | Let's make a generic template subroutine using this three-point central difference | |
238 | rule. In particular: | |
239 | */ | |
240 | ||
241 | template<typename value_type, typename function_type> | |
242 | value_type derivative(const value_type x, const value_type dx, function_type func) | |
243 | { | |
244 | // Compute d/dx[func(*first)] using a three-point | |
245 | // central difference rule of O(dx^6). | |
246 | ||
247 | const value_type dx1 = dx; | |
248 | const value_type dx2 = dx1 * 2; | |
249 | const value_type dx3 = dx1 * 3; | |
250 | ||
251 | const value_type m1 = (func(x + dx1) - func(x - dx1)) / 2; | |
252 | const value_type m2 = (func(x + dx2) - func(x - dx2)) / 4; | |
253 | const value_type m3 = (func(x + dx3) - func(x - dx3)) / 6; | |
254 | ||
255 | const value_type fifteen_m1 = 15 * m1; | |
256 | const value_type six_m2 = 6 * m2; | |
257 | const value_type ten_dx1 = 10 * dx1; | |
258 | ||
259 | return ((fifteen_m1 - six_m2) + m3) / ten_dx1; | |
260 | } | |
261 | ||
262 | /*`The `derivative()` template function can be used to compute the first derivative | |
263 | of any function to ['O(dx[super 6])]. For example, consider the first derivative of ['sin(x)] evaluated | |
264 | at ['x = [pi]/3]. In other words, | |
265 | ||
266 | [equation floating_point_eg3] | |
267 | ||
268 | The code below computes the derivative in Equation 3 for float, double and boost's | |
269 | multiple-precision type cpp_dec_float_50. | |
270 | */ | |
271 | ||
272 | //] | |
273 | ||
274 | //[GI1 | |
275 | ||
276 | /*` | |
277 | Similar to the generic derivative example, we can calculate integrals in a similar manner: | |
278 | */ | |
279 | ||
280 | template<typename value_type, typename function_type> | |
281 | inline value_type integral(const value_type a, | |
282 | const value_type b, | |
283 | const value_type tol, | |
284 | function_type func) | |
285 | { | |
286 | unsigned n = 1U; | |
287 | ||
288 | value_type h = (b - a); | |
289 | value_type I = (func(a) + func(b)) * (h / 2); | |
290 | ||
291 | for(unsigned k = 0U; k < 8U; k++) | |
292 | { | |
293 | h /= 2; | |
294 | ||
295 | value_type sum(0); | |
296 | for(unsigned j = 1U; j <= n; j++) | |
297 | { | |
298 | sum += func(a + (value_type((j * 2) - 1) * h)); | |
299 | } | |
300 | ||
301 | const value_type I0 = I; | |
302 | I = (I / 2) + (h * sum); | |
303 | ||
304 | const value_type ratio = I0 / I; | |
305 | const value_type delta = ratio - 1; | |
306 | const value_type delta_abs = ((delta < 0) ? -delta : delta); | |
307 | ||
308 | if((k > 1U) && (delta_abs < tol)) | |
309 | { | |
310 | break; | |
311 | } | |
312 | ||
313 | n *= 2U; | |
314 | } | |
315 | ||
316 | return I; | |
317 | } | |
318 | ||
319 | /*` | |
320 | The following sample program shows how the function can be called, we begin | |
321 | by defining a function object, which when integrated should yield the Bessel J | |
322 | function: | |
323 | */ | |
324 | ||
325 | template<typename value_type> | |
326 | class cyl_bessel_j_integral_rep | |
327 | { | |
328 | public: | |
329 | cyl_bessel_j_integral_rep(const unsigned N, | |
330 | const value_type& X) : n(N), x(X) { } | |
331 | ||
332 | value_type operator()(const value_type& t) const | |
333 | { | |
334 | // pi * Jn(x) = Int_0^pi [cos(x * sin(t) - n*t) dt] | |
335 | return cos(x * sin(t) - (n * t)); | |
336 | } | |
337 | ||
338 | private: | |
339 | const unsigned n; | |
340 | const value_type x; | |
341 | }; | |
342 | ||
343 | ||
344 | //] | |
345 | ||
346 | //[POLY | |
347 | ||
348 | /*` | |
349 | In this example we'll look at polynomial evaluation, this is not only an important | |
350 | use case, but it's one that `number` performs particularly well at because the | |
351 | expression templates ['completely eliminate all temporaries] from a | |
352 | [@http://en.wikipedia.org/wiki/Horner%27s_method Horner polynomial | |
353 | evaluation scheme]. | |
354 | ||
355 | The following code evaluates `sin(x)` as a polynomial, accurate to at least 64 decimal places: | |
356 | ||
357 | */ | |
358 | ||
359 | using boost::multiprecision::cpp_dec_float; | |
360 | typedef boost::multiprecision::number<cpp_dec_float<64> > mp_type; | |
361 | ||
362 | mp_type mysin(const mp_type& x) | |
363 | { | |
364 | // Approximation of sin(x * pi/2) for -1 <= x <= 1, using an order 63 polynomial. | |
365 | static const std::array<mp_type, 32U> coefs = | |
366 | {{ | |
367 | mp_type("+1.5707963267948966192313216916397514420985846996875529104874722961539082031431044993140174126711"), //"), | |
368 | mp_type("-0.64596409750624625365575656389794573337969351178927307696134454382929989411386887578263960484"), // ^3 | |
369 | mp_type("+0.07969262624616704512050554949047802252091164235106119545663865720995702920146198554317279"), // ^5 | |
370 | mp_type("-0.0046817541353186881006854639339534378594950280185010575749538605102665157913157426229824"), // ^7 | |
371 | mp_type("+0.00016044118478735982187266087016347332970280754062061156858775174056686380286868007443"), // ^9 | |
372 | mp_type("-3.598843235212085340458540018208389404888495232432127661083907575106196374913134E-6"), // ^11 | |
373 | mp_type("+5.692172921967926811775255303592184372902829756054598109818158853197797542565E-8"), // ^13 | |
374 | mp_type("-6.688035109811467232478226335783138689956270985704278659373558497256423498E-10"), // ^15 | |
375 | mp_type("+6.066935731106195667101445665327140070166203261129845646380005577490472E-12"), // ^17 | |
376 | mp_type("-4.377065467313742277184271313776319094862897030084226361576452003432E-14"), // ^19 | |
377 | mp_type("+2.571422892860473866153865950420487369167895373255729246889168337E-16"), // ^21 | |
378 | mp_type("-1.253899540535457665340073300390626396596970180355253776711660E-18"), // ^23 | |
379 | mp_type("+5.15645517658028233395375998562329055050964428219501277474E-21"), // ^25 | |
380 | mp_type("-1.812399312848887477410034071087545686586497030654642705E-23"), // ^27 | |
381 | mp_type("+5.50728578652238583570585513920522536675023562254864E-26"), // ^29 | |
382 | mp_type("-1.461148710664467988723468673933026649943084902958E-28"), // ^31 | |
383 | mp_type("+3.41405297003316172502972039913417222912445427E-31"), // ^33 | |
384 | mp_type("-7.07885550810745570069916712806856538290251E-34"), // ^35 | |
385 | mp_type("+1.31128947968267628970845439024155655665E-36"), // ^37 | |
386 | mp_type("-2.18318293181145698535113946654065918E-39"), // ^39 | |
387 | mp_type("+3.28462680978498856345937578502923E-42"), // ^41 | |
388 | mp_type("-4.48753699028101089490067137298E-45"), // ^43 | |
389 | mp_type("+5.59219884208696457859353716E-48"), // ^45 | |
390 | mp_type("-6.38214503973500471720565E-51"), // ^47 | |
391 | mp_type("+6.69528558381794452556E-54"), // ^49 | |
392 | mp_type("-6.47841373182350206E-57"), // ^51 | |
393 | mp_type("+5.800016389666445E-60"), // ^53 | |
394 | mp_type("-4.818507347289E-63"), // ^55 | |
395 | mp_type("+3.724683686E-66"), // ^57 | |
396 | mp_type("-2.6856479E-69"), // ^59 | |
397 | mp_type("+1.81046E-72"), // ^61 | |
398 | mp_type("-1.133E-75"), // ^63 | |
399 | }}; | |
400 | ||
401 | const mp_type v = x * 2 / boost::math::constants::pi<mp_type>(); | |
402 | const mp_type x2 = (v * v); | |
403 | // | |
404 | // Polynomial evaluation follows, if mp_type allocates memory then | |
405 | // just one such allocation occurs - to initialize the variable "sum" - | |
406 | // and no temporaries are created at all. | |
407 | // | |
408 | const mp_type sum = ((((((((((((((((((((((((((((((( + coefs[31U] | |
409 | * x2 + coefs[30U]) | |
410 | * x2 + coefs[29U]) | |
411 | * x2 + coefs[28U]) | |
412 | * x2 + coefs[27U]) | |
413 | * x2 + coefs[26U]) | |
414 | * x2 + coefs[25U]) | |
415 | * x2 + coefs[24U]) | |
416 | * x2 + coefs[23U]) | |
417 | * x2 + coefs[22U]) | |
418 | * x2 + coefs[21U]) | |
419 | * x2 + coefs[20U]) | |
420 | * x2 + coefs[19U]) | |
421 | * x2 + coefs[18U]) | |
422 | * x2 + coefs[17U]) | |
423 | * x2 + coefs[16U]) | |
424 | * x2 + coefs[15U]) | |
425 | * x2 + coefs[14U]) | |
426 | * x2 + coefs[13U]) | |
427 | * x2 + coefs[12U]) | |
428 | * x2 + coefs[11U]) | |
429 | * x2 + coefs[10U]) | |
430 | * x2 + coefs[9U]) | |
431 | * x2 + coefs[8U]) | |
432 | * x2 + coefs[7U]) | |
433 | * x2 + coefs[6U]) | |
434 | * x2 + coefs[5U]) | |
435 | * x2 + coefs[4U]) | |
436 | * x2 + coefs[3U]) | |
437 | * x2 + coefs[2U]) | |
438 | * x2 + coefs[1U]) | |
439 | * x2 + coefs[0U]) | |
440 | * v; | |
441 | ||
442 | return sum; | |
443 | } | |
444 | ||
445 | /*` | |
446 | Calling the function like so: | |
447 | ||
448 | mp_type pid4 = boost::math::constants::pi<mp_type>() / 4; | |
449 | std::cout << std::setprecision(std::numeric_limits< ::mp_type>::digits10) << std::scientific; | |
450 | std::cout << mysin(pid4) << std::endl; | |
451 | ||
452 | Yields the expected output: | |
453 | ||
454 | [pre 7.0710678118654752440084436210484903928483593768847403658833986900e-01] | |
455 | ||
456 | */ | |
457 | ||
458 | //] | |
459 | ||
460 | ||
461 | int main() | |
462 | { | |
463 | using namespace boost::multiprecision; | |
464 | std::cout << std::scientific << std::setprecision(std::numeric_limits<double>::digits10); | |
465 | std::cout << JEL1(2.5, 0.5) << std::endl; | |
466 | std::cout << std::scientific << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10); | |
467 | std::cout << JEL2(cpp_dec_float_50(2.5), cpp_dec_float_50(0.5)) << std::endl; | |
468 | std::cout << std::scientific << std::setprecision(std::numeric_limits<double>::digits10); | |
469 | std::cout << JEL3(2.5, 0.5) << std::endl; | |
470 | std::cout << std::scientific << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10); | |
471 | std::cout << JEL3(cpp_dec_float_50(2.5), cpp_dec_float_50(0.5)) << std::endl; | |
472 | std::cout << std::scientific << std::setprecision(std::numeric_limits<cpp_dec_float_100>::digits10); | |
473 | std::cout << JEL4(cpp_dec_float_100(2) + 0.5, cpp_dec_float_100(0.5)) << std::endl; | |
474 | ||
475 | //[AOS2 | |
476 | ||
477 | /*=#include <iostream> | |
478 | #include <iomanip> | |
479 | #include <boost/multiprecision/cpp_dec_float.hpp> | |
480 | ||
481 | using boost::multiprecision::cpp_dec_float_50; | |
482 | ||
483 | int main(int, char**) | |
484 | {*/ | |
485 | const float r_f(float(123) / 100); | |
486 | const float a_f = area_of_a_circle(r_f); | |
487 | ||
488 | const double r_d(double(123) / 100); | |
489 | const double a_d = area_of_a_circle(r_d); | |
490 | ||
491 | const cpp_dec_float_50 r_mp(cpp_dec_float_50(123) / 100); | |
492 | const cpp_dec_float_50 a_mp = area_of_a_circle(r_mp); | |
493 | ||
494 | // 4.75292 | |
495 | std::cout | |
496 | << std::setprecision(std::numeric_limits<float>::digits10) | |
497 | << a_f | |
498 | << std::endl; | |
499 | ||
500 | // 4.752915525616 | |
501 | std::cout | |
502 | << std::setprecision(std::numeric_limits<double>::digits10) | |
503 | << a_d | |
504 | << std::endl; | |
505 | ||
506 | // 4.7529155256159981904701331745635599135018975843146 | |
507 | std::cout | |
508 | << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10) | |
509 | << a_mp | |
510 | << std::endl; | |
511 | /*=}*/ | |
512 | ||
513 | //] | |
514 | ||
515 | //[ND2 | |
516 | /*= | |
517 | #include <iostream> | |
518 | #include <iomanip> | |
519 | #include <boost/multiprecision/cpp_dec_float.hpp> | |
520 | #include <boost/math/constants/constants.hpp> | |
521 | ||
522 | ||
523 | int main(int, char**) | |
524 | {*/ | |
525 | using boost::math::constants::pi; | |
526 | using boost::multiprecision::cpp_dec_float_50; | |
527 | // | |
528 | // We'll pass a function pointer for the function object passed to derivative, | |
529 | // the typecast is needed to select the correct overload of std::sin: | |
530 | // | |
531 | const float d_f = derivative( | |
532 | pi<float>() / 3, | |
533 | 0.01F, | |
534 | static_cast<float(*)(float)>(std::sin) | |
535 | ); | |
536 | ||
537 | const double d_d = derivative( | |
538 | pi<double>() / 3, | |
539 | 0.001, | |
540 | static_cast<double(*)(double)>(std::sin) | |
541 | ); | |
542 | // | |
543 | // In the cpp_dec_float_50 case, the sin function is multiply overloaded | |
544 | // to handle expression templates etc. As a result it's hard to take its | |
545 | // address without knowing about its implementation details. We'll use a | |
546 | // C++11 lambda expression to capture the call. | |
547 | // We also need a typecast on the first argument so we don't accidentally pass | |
548 | // an expression template to a template function: | |
549 | // | |
550 | const cpp_dec_float_50 d_mp = derivative( | |
551 | cpp_dec_float_50(pi<cpp_dec_float_50>() / 3), | |
552 | cpp_dec_float_50(1.0E-9), | |
553 | [](const cpp_dec_float_50& x) -> cpp_dec_float_50 | |
554 | { | |
555 | return sin(x); | |
556 | } | |
557 | ); | |
558 | ||
559 | // 5.000029e-001 | |
560 | std::cout | |
561 | << std::setprecision(std::numeric_limits<float>::digits10) | |
562 | << d_f | |
563 | << std::endl; | |
564 | ||
565 | // 4.999999999998876e-001 | |
566 | std::cout | |
567 | << std::setprecision(std::numeric_limits<double>::digits10) | |
568 | << d_d | |
569 | << std::endl; | |
570 | ||
571 | // 4.99999999999999999999999999999999999999999999999999e-01 | |
572 | std::cout | |
573 | << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10) | |
574 | << d_mp | |
575 | << std::endl; | |
576 | //=} | |
577 | ||
578 | /*` | |
579 | The expected value of the derivative is 0.5. This central difference rule in this | |
580 | example is ill-conditioned, meaning it suffers from slight loss of precision. With that | |
581 | in mind, the results agree with the expected value of 0.5.*/ | |
582 | ||
583 | //] | |
584 | ||
585 | //[ND3 | |
586 | ||
587 | /*` | |
588 | We can take this a step further and use our derivative function to compute | |
589 | a partial derivative. For example if we take the incomplete gamma function | |
590 | ['P(a, z)], and take the derivative with respect to /z/ at /(2,2)/ then we | |
591 | can calculate the result as shown below, for good measure we'll compare with | |
592 | the "correct" result obtained from a call to ['gamma_p_derivative], the results | |
593 | agree to approximately 44 digits: | |
594 | */ | |
595 | ||
596 | cpp_dec_float_50 gd = derivative( | |
597 | cpp_dec_float_50(2), | |
598 | cpp_dec_float_50(1.0E-9), | |
599 | [](const cpp_dec_float_50& x) ->cpp_dec_float_50 | |
600 | { | |
601 | return boost::math::gamma_p(2, x); | |
602 | } | |
603 | ); | |
604 | // 2.70670566473225383787998989944968806815263091819151e-01 | |
605 | std::cout | |
606 | << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10) | |
607 | << gd | |
608 | << std::endl; | |
609 | // 2.70670566473225383787998989944968806815253190143120e-01 | |
610 | std::cout << boost::math::gamma_p_derivative(cpp_dec_float_50(2), cpp_dec_float_50(2)) << std::endl; | |
611 | //] | |
612 | ||
613 | //[GI2 | |
614 | ||
615 | /* The function can now be called as follows: */ | |
616 | /*=int main(int, char**) | |
617 | {*/ | |
618 | using boost::math::constants::pi; | |
619 | typedef boost::multiprecision::cpp_dec_float_50 mp_type; | |
620 | ||
621 | const float j2_f = | |
622 | integral(0.0F, | |
623 | pi<float>(), | |
624 | 0.01F, | |
625 | cyl_bessel_j_integral_rep<float>(2U, 1.23F)) / pi<float>(); | |
626 | ||
627 | const double j2_d = | |
628 | integral(0.0, | |
629 | pi<double>(), | |
630 | 0.0001, | |
631 | cyl_bessel_j_integral_rep<double>(2U, 1.23)) / pi<double>(); | |
632 | ||
633 | const mp_type j2_mp = | |
634 | integral(mp_type(0), | |
635 | pi<mp_type>(), | |
636 | mp_type(1.0E-20), | |
637 | cyl_bessel_j_integral_rep<mp_type>(2U, mp_type(123) / 100)) / pi<mp_type>(); | |
638 | ||
639 | // 0.166369 | |
640 | std::cout | |
641 | << std::setprecision(std::numeric_limits<float>::digits10) | |
642 | << j2_f | |
643 | << std::endl; | |
644 | ||
645 | // 0.166369383786814 | |
646 | std::cout | |
647 | << std::setprecision(std::numeric_limits<double>::digits10) | |
648 | << j2_d | |
649 | << std::endl; | |
650 | ||
651 | // 0.16636938378681407351267852431513159437103348245333 | |
652 | std::cout | |
653 | << std::setprecision(std::numeric_limits<mp_type>::digits10) | |
654 | << j2_mp | |
655 | << std::endl; | |
656 | ||
657 | // | |
658 | // Print true value for comparison: | |
659 | // 0.166369383786814073512678524315131594371033482453329 | |
660 | std::cout << boost::math::cyl_bessel_j(2, mp_type(123) / 100) << std::endl; | |
661 | //=} | |
662 | ||
663 | //] | |
664 | ||
665 | std::cout << std::setprecision(std::numeric_limits< ::mp_type>::digits10) << std::scientific; | |
666 | std::cout << mysin(boost::math::constants::pi< ::mp_type>() / 4) << std::endl; | |
667 | std::cout << boost::multiprecision::sin(boost::math::constants::pi< ::mp_type>() / 4) << std::endl; | |
668 | ||
669 | return 0; | |
670 | } | |
671 | ||
672 | /* | |
673 | ||
674 | Program output: | |
675 | ||
676 | 9.822663964796047e-001 | |
677 | 9.82266396479604757017335009796882833995903762577173e-01 | |
678 | 9.822663964796047e-001 | |
679 | 9.82266396479604757017335009796882833995903762577173e-01 | |
680 | 9.8226639647960475701733500979688283399590376257717309069410413822165082248153638454147004236848917775e-01 | |
681 | 4.752916e+000 | |
682 | 4.752915525615998e+000 | |
683 | 4.75291552561599819047013317456355991350189758431460e+00 | |
684 | 5.000029e-001 | |
685 | 4.999999999998876e-001 | |
686 | 4.99999999999999999999999999999999999999999999999999e-01 | |
687 | 2.70670566473225383787998989944968806815263091819151e-01 | |
688 | 2.70670566473225383787998989944968806815253190143120e-01 | |
689 | 7.0710678118654752440084436210484903928483593768847403658833986900e-01 | |
690 | 7.0710678118654752440084436210484903928483593768847403658833986900e-01 | |
691 | */ | |
692 | ||
693 | #else | |
694 | ||
695 | int main() { return 0; } | |
696 | ||
697 | #endif |