2

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.

CuriousQ
  • 31
  • 5
  • what do you mean by "it gets stuck"? is there an error or exception or does it just hang? – jwezorek Aug 04 '21 at 18:18
  • @jwezorek thanks for your reply. It compiles without any warnings. It hangs when calculating the exponent. It does something, but it has never computed the exponent. Even when setting the number of digits to 16 (but still using mpc_complex_backened). – CuriousQ Aug 04 '21 at 19:15
  • can you step into it in a debugger? – jwezorek Aug 04 '21 at 19:18
  • Yes. I tried to go ``into'' as far as it is possible. Debugger threw an error "unable to resolve non-existing file.../src/init2.c". It threw this error when calling a function mpc_init2 from mpc_complex_imp(). It can be stepped over. When stepping over, it hangs withing a template BOOST_MP_CXX14_CONSTEXPR number(...). Maybe there is a better way to collect all this sort of information? I can try to find out. – CuriousQ Aug 04 '21 at 19:34
  • intereseting that it is having problems in an `init` function. do you know if you have GMP installed? I'm not on Linux so cant test this but visiual studio fails to build your code for me because I do not have GMP. Perhaps you dont either but it actually compiles with your set up? – jwezorek Aug 04 '21 at 19:41
  • I think I have installed gmp. At least I have a header file gmp.h. Also, on the start of the debugger, it loads libgmp.so.10. So i think it is ananother issue. Let me try to collect call stack. – CuriousQ Aug 04 '21 at 19:50
  • So it is something to do with libgmp. Namely, it hangs on calling ___gmpn_sub_n_coreisbr(). I can try to follow your suggestion,and check if libgmp.so.10 exists somewhere on my system. – CuriousQ Aug 04 '21 at 20:02
  • yeah it has something to do with that dependency. I can't help you more as this seems like a Linux packages issue so I don't know – jwezorek Aug 04 '21 at 20:06
  • @jwezorek thanks anyway! One thing though. I tried to use cpp_complex_backened, and the code worked just fine. Though according to Boost webpage, it much slower than mpc... – CuriousQ Aug 04 '21 at 20:16

1 Answers1

1

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.

CuriousQ
  • 31
  • 5