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_
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>
13 #if !defined(BOOST_NO_CXX11_HDR_ARRAY) && !defined(BOOST_NO_CXX11_LAMBDAS)
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
22 [:['a = [pi] * r[super 2]]]
24 The area of a circle can be computed in generic programming using Boost.Math
25 for the constant [pi] as shown below:
29 //=#include <boost/math/constants/constants.hpp>
32 inline T
area_of_a_circle(T r
)
34 using boost::math::constants::pi
;
35 return pi
<T
>() * r
* r
;
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
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.
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.*/
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.
66 The Jahnke-Emden Lambda function is defined by the equation:
68 [:['JahnkeEmden(v, z) = [Gamma](v+1) * J[sub v](z) / (z / 2)[super v]]]
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:
75 double JEL1(double v
, double z
)
77 return boost::math::tgamma(v
+ 1) * boost::math::cyl_bessel_j(v
, z
) / std::pow(z
/ 2, v
);
81 Calling this function as:
83 std::cout << std::scientific << std::setprecision(std::numeric_limits<double>::digits10);
84 std::cout << JEL1(2.5, 0.5) << std::endl;
88 [pre 9.822663964796047e-001]
90 Now let's implement the function again, but this time using the multiprecision type
91 `cpp_dec_float_50` as the argument type:
95 boost::multiprecision::cpp_dec_float_50
96 JEL2(boost::multiprecision::cpp_dec_float_50 v
, boost::multiprecision::cpp_dec_float_50 z
)
98 return boost::math::tgamma(v
+ 1) * boost::math::cyl_bessel_j(v
, z
) / boost::multiprecision::pow(z
/ 2, v
);
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.
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
112 Here's an example of how the function may be called:
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;
119 [pre 9.82266396479604757017335009796882833995903762577173e-01]
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:
126 template <class Float
>
127 Float
JEL3(Float v
, Float z
)
130 return boost::math::tgamma(v
+ 1) * boost::math::cyl_bessel_j(v
, z
) / pow(z
/ 2, v
);
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.
142 Let's call the function with both `double` and multiprecision arguments:
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;
152 9.822663964796047e-001
153 9.82266396479604757017335009796882833995903762577173e-01
156 Unfortunately there is a problem with this version: if we were to call it like this:
158 boost::multiprecision::cpp_dec_float_50 v(2), z(0.5);
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
165 JEL(cpp_dec_float_50(v + 0.5), z);
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:
173 template <class Float1
, class Float2
>
174 typename
boost::math::tools::promote_args
<Float1
, Float2
>::type
175 JEL4(Float1 v
, Float2 z
)
178 return boost::math::tgamma(v
+ 1) * boost::math::cyl_bessel_j(v
, z
) / pow(z
/ 2, v
);
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.
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;
191 And get 100 digits of output:
193 [pre 9.8226639647960475701733500979688283399590376257717309069410413822165082248153638454147004236848917775e-01]
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.
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:
202 template <class Float1, class Float2>
203 typename boost::math::tools::promote_args<Float1, Float2>::type
204 JEL5(Float1 v, Float2 z)
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;
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]]:
224 [equation floating_point_eg1]
226 Where the difference terms ['m[sub n]] are given by:
228 [equation floating_point_eg2]
230 and ['dx] is the step-size of the derivative.
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.
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
241 template<typename value_type
, typename function_type
>
242 value_type
derivative(const value_type x
, const value_type dx
, function_type func
)
244 // Compute d/dx[func(*first)] using a three-point
245 // central difference rule of O(dx^6).
247 const value_type dx1
= dx
;
248 const value_type dx2
= dx1
* 2;
249 const value_type dx3
= dx1
* 3;
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;
255 const value_type fifteen_m1
= 15 * m1
;
256 const value_type six_m2
= 6 * m2
;
257 const value_type ten_dx1
= 10 * dx1
;
259 return ((fifteen_m1
- six_m2
) + m3
) / ten_dx1
;
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,
266 [equation floating_point_eg3]
268 The code below computes the derivative in Equation 3 for float, double and boost's
269 multiple-precision type cpp_dec_float_50.
277 Similar to the generic derivative example, we can calculate integrals in a similar manner:
280 template<typename value_type
, typename function_type
>
281 inline value_type
integral(const value_type a
,
283 const value_type tol
,
288 value_type h
= (b
- a
);
289 value_type I
= (func(a
) + func(b
)) * (h
/ 2);
291 for(unsigned k
= 0U; k
< 8U; k
++)
296 for(unsigned j
= 1U; j
<= n
; j
++)
298 sum
+= func(a
+ (value_type((j
* 2) - 1) * h
));
301 const value_type I0
= I
;
302 I
= (I
/ 2) + (h
* sum
);
304 const value_type ratio
= I0
/ I
;
305 const value_type delta
= ratio
- 1;
306 const value_type delta_abs
= ((delta
< 0) ? -delta
: delta
);
308 if((k
> 1U) && (delta_abs
< tol
))
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
325 template<typename value_type
>
326 class cyl_bessel_j_integral_rep
329 cyl_bessel_j_integral_rep(const unsigned N
,
330 const value_type
& X
) : n(N
), x(X
) { }
332 value_type
operator()(const value_type
& t
) const
334 // pi * Jn(x) = Int_0^pi [cos(x * sin(t) - n*t) dt]
335 return cos(x
* sin(t
) - (n
* t
));
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
355 The following code evaluates `sin(x)` as a polynomial, accurate to at least 64 decimal places:
359 using boost::multiprecision::cpp_dec_float
;
360 typedef boost::multiprecision::number
<cpp_dec_float
<64> > mp_type
;
362 mp_type
mysin(const mp_type
& x
)
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
=
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
401 const mp_type v
= x
* 2 / boost::math::constants::pi
<mp_type
>();
402 const mp_type x2
= (v
* v
);
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.
408 const mp_type sum
= ((((((((((((((((((((((((((((((( + coefs
[31U]
446 Calling the function like so:
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;
452 Yields the expected output:
454 [pre 7.0710678118654752440084436210484903928483593768847403658833986900e-01]
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
;
477 /*=#include <iostream>
479 #include <boost/multiprecision/cpp_dec_float.hpp>
481 using boost::multiprecision::cpp_dec_float_50;
483 int main(int, char**)
485 const float r_f(float(123) / 100);
486 const float a_f
= area_of_a_circle(r_f
);
488 const double r_d(double(123) / 100);
489 const double a_d
= area_of_a_circle(r_d
);
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
);
496 << std::setprecision(std::numeric_limits
<float>::digits10
)
502 << std::setprecision(std::numeric_limits
<double>::digits10
)
506 // 4.7529155256159981904701331745635599135018975843146
508 << std::setprecision(std::numeric_limits
<cpp_dec_float_50
>::digits10
)
519 #include <boost/multiprecision/cpp_dec_float.hpp>
520 #include <boost/math/constants/constants.hpp>
523 int main(int, char**)
525 using boost::math::constants::pi
;
526 using boost::multiprecision::cpp_dec_float_50
;
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:
531 const float d_f
= derivative(
534 static_cast<float(*)(float)>(std::sin
)
537 const double d_d
= derivative(
540 static_cast<double(*)(double)>(std::sin
)
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:
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
561 << std::setprecision(std::numeric_limits
<float>::digits10
)
565 // 4.999999999998876e-001
567 << std::setprecision(std::numeric_limits
<double>::digits10
)
571 // 4.99999999999999999999999999999999999999999999999999e-01
573 << std::setprecision(std::numeric_limits
<cpp_dec_float_50
>::digits10
)
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.*/
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:
596 cpp_dec_float_50 gd
= derivative(
598 cpp_dec_float_50(1.0E-9),
599 [](const cpp_dec_float_50
& x
) ->cpp_dec_float_50
601 return boost::math::gamma_p(2, x
);
604 // 2.70670566473225383787998989944968806815263091819151e-01
606 << std::setprecision(std::numeric_limits
<cpp_dec_float_50
>::digits10
)
609 // 2.70670566473225383787998989944968806815253190143120e-01
610 std::cout
<< boost::math::gamma_p_derivative(cpp_dec_float_50(2), cpp_dec_float_50(2)) << std::endl
;
615 /* The function can now be called as follows: */
616 /*=int main(int, char**)
618 using boost::math::constants::pi
;
619 typedef boost::multiprecision::cpp_dec_float_50 mp_type
;
625 cyl_bessel_j_integral_rep
<float>(2U, 1.23F
)) / pi
<float>();
631 cyl_bessel_j_integral_rep
<double>(2U, 1.23)) / pi
<double>();
633 const mp_type j2_mp
=
637 cyl_bessel_j_integral_rep
<mp_type
>(2U, mp_type(123) / 100)) / pi
<mp_type
>();
641 << std::setprecision(std::numeric_limits
<float>::digits10
)
647 << std::setprecision(std::numeric_limits
<double>::digits10
)
651 // 0.16636938378681407351267852431513159437103348245333
653 << std::setprecision(std::numeric_limits
<mp_type
>::digits10
)
658 // Print true value for comparison:
659 // 0.166369383786814073512678524315131594371033482453329
660 std::cout
<< boost::math::cyl_bessel_j(2, mp_type(123) / 100) << std::endl
;
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
;
676 9.822663964796047e-001
677 9.82266396479604757017335009796882833995903762577173e-01
678 9.822663964796047e-001
679 9.82266396479604757017335009796882833995903762577173e-01
680 9.8226639647960475701733500979688283399590376257717309069410413822165082248153638454147004236848917775e-01
682 4.752915525615998e+000
683 4.75291552561599819047013317456355991350189758431460e+00
685 4.999999999998876e-001
686 4.99999999999999999999999999999999999999999999999999e-01
687 2.70670566473225383787998989944968806815263091819151e-01
688 2.70670566473225383787998989944968806815253190143120e-01
689 7.0710678118654752440084436210484903928483593768847403658833986900e-01
690 7.0710678118654752440084436210484903928483593768847403658833986900e-01
695 int main() { return 0; }