为什么 boost::multiprecision::exp 在评估复值积分时卡住？

Question

我是 Boost C++ 库的新手，在使用它们时自然遇到了很多问题（由于缺乏知识和示例:)

其中一个问题来自以下一段代码


    #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;
    }

当代码到达 mpc_type ex2 = mpns::exp(-z*x2) ; 时卡住了。即，它在计算指数时挂起。我花了一些时间试图找出问题所在，但找不到解决方案。

我做了几个测试。例如，<boost/multiprecision/mpfr.hpp> 工作得很好。，即被积函数的实值版本可以积分到任意精度。我测试了代码的 mpfr 类型版本，直到 boost::multiprecision::mpfr_float_backend <2000>.

调用 lambda 函数 f 是可能的，它 returns 正确的数字。

我使用了不同的被积函数，例如z*x 、z*tgamma(x) 并且程序运行良好，z 和 x 的定义相同，您可以在上面的示例中找到。

我使用 Tumbleweed 提供的最新版本的 Boost 库，即 boost_1_76_0-gnu-mpich-hpc-devel

编译器：g++

cpp标准：gnu++17

这个问题的根源是什么？

抱歉问题很长。

Answer 1

感谢来自 Boost 的 John Maddock，问题得以解决。也就是说，约翰指出，尽管对指数的值施加了限制，但结果变得非常小。当mpc_exp接近这样一个极小的值时，它会越来越慢。

发生这种情况的原因在 https://www.boost.org/doc/libs/1_76_0/libs/math/doc/html/math_toolkit/double_exponential/de_caveats.html

中有所描述

此问题的解决方法是使用不同的约束。我按照约翰的建议使用了以下内容（使用与以前相同的符号）


        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)) ;
        }

在选定的指数范围内，原则上可以积分到 3000 位。

为什么 boost::multiprecision::exp 在评估复值积分时卡住？

Why boost::multiprecision::exp gets stuck when evaluating complex-valued integral?

c++

boost

boost-multiprecision