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>
area_of_a_circle(T r)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 
JEL1(double v,double z)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
JEL2(boost::multiprecision::cpp_dec_float_50 v,boost::multiprecision::cpp_dec_float_50 z)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>
JEL3(Float v,Float z)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
JEL4(Float1 v,Float2 z)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>
integral(const value_type a,const value_type b,const value_type tol,function_type func)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:
cyl_bessel_j_integral_rep(const unsigned N,const value_type & X)329    cyl_bessel_j_integral_rep(const unsigned N,
330       const value_type& X) : n(N), x(X) { }
331 
operator ()(const value_type & t) const332    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 
mysin(const mp_type & x)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 
main()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 
main()695 int main() { return 0; }
696 
697 #endif
698