I am new to Boost C++ Libraries, and naturally I have encountered many problems when using them (due to lack of knowledge and examples available :)
One of these problems comes from the following piece of code
#include <iostream>
#include <boost/math/constants/constants.hpp>
#include <boost/math/quadrature/exp_sinh.hpp>
#include <boost/multiprecision/mpc.hpp>
#include <boost/multiprecision/mpfr.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
namespace mpns = boost::multiprecision;
typedef mpns::number<boost::multiprecision::mpc_complex_backend <50> > mpc_type ;
typedef mpc_type::value_type mp_type ;
int main()
{
using boost::math::constants::root_pi ;
mpc_type z{mp_type(1),mp_type(1)} ;
auto f = [&](mp_type x)
{
//actually I did not test if all these conditions are needed...
if (boost::math::fpclassify<mp_type> (x) == FP_ZERO)
{
return mpc_type(mp_type(1)) ;
};
if (boost::math::fpclassify<mp_type> (x) == FP_INFINITE)
{
return mpc_type(mp_type(0)) ;
};
mp_type x2 = mpns::pow(x,2U) ;
if (boost::math::fpclassify<mp_type> (x2) == FP_ZERO)
{
return mpc_type(mp_type(1)) ;
};
if (boost::math::fpclassify<mp_type> (x2) == FP_INFINITE)
{
return mpc_type(mp_type(0)) ;
};
mpc_type ex2 = mpns::exp(-z*x2) ;
if (boost::math::fpclassify<mpc_type> (ex2) == FP_ZERO)
{
return mpc_type(mp_type(0)) ;
};
return ex2 ;
} ;
mp_type termination = boost::math::tools::root_epsilon <mp_type> () ;
mp_type error;
mp_type L1;
size_t max_halvings = 12;
boost::math::quadrature::exp_sinh<mp_type> integrator(max_halvings) ;
mpc_type res = integrator.integrate(f,termination,&error,&L1) ;
mpc_type div = 2U*mpns::sqrt(z) ;
mpc_type result = (mpc_type(root_pi<mp_type> ())/div) - res ;
return 0;
}
When the code reaches mpc_type ex2 = mpns::exp(-z*x2) ;it gets stuck. Namely, it hangs when calculating the exponent. I spent some time trying to find out what is wrong but could not find a solution.
I did several tests.
For instance, <boost/multiprecision/mpfr.hpp> works just fine.,i.e. a real-valued version of the integrand can be integrated to arbitrary precision. I tested the mpfr-type version of the code up to boost::multiprecision::mpfr_float_backend <2000>.
Calling the lambda-function f is possible and it returns correct numbers.
I used different integrands, e.g. z*x , z*tgamma(x) and the program worked fine, with the same definitions of z and x which you can find in the example above.
I use the latest version of the Boost libraries as provided by Tumbleweed, i.e. boost_1_76_0-gnu-mpich-hpc-devel
Compiler : g++
cppStandard : gnu++17
What could be a source of this problem?
Sorry for a lengthy question.
Thanks to John Maddock from Boost, it was possible to resolve the problem. Namely, John pointed out that, despite constraints imposed on the value of the exponent, the result was getting extremely small. When mpc_exp is approaching such an extremely small value, it gets slower and slower.
The reason for this to be happening is described in, for instance, https://www.boost.org/doc/libs/1_76_0/libs/math/doc/html/math_toolkit/double_exponential/de_caveats.html
A workaround for this problem is to use different constraints. I used the following, as suggested by John (using the same notation as before)
mp_type x2 = mpns::pow(x,2U) ;
int minexpx2 = -10000 ;
int maxexpx2 = 10000;
int exponentx2 ;
mp_type tmp = frexp(x2,&exponentx2) ;
if (exponentx2 <= minexpx2)
{
return mpc_type(mp_type(1)) ;
}
else if (exponentx2 >= maxexpx2)
{
return mpc_type(mp_type(0)) ;
}
With the chosen range of exponents, it is possible to perform integration, in principle, to 3000 digits.