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

Posted 2023-02-21

【中文标题】为啥在评估复值积分时 boost::multiprecision::exp 会卡住？

【英文标题】：Why boost::multiprecision::exp gets stuck when evaluating complex-valued integral?为什么在评估复值积分时 boost::multiprecision::exp 会卡住？

【发布时间】：2021-08-04 13:23:11

【问题描述】：

我是 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 zmp_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-type 版本，直到boost::multiprecision::mpfr_float_backend <2000>。

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

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

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

编译器：g++

cpp标准：gnu++17

这个问题的根源是什么？

抱歉，问了一个冗长的问题。

【问题讨论】：

“卡住”是什么意思？是有错误或异常还是只是挂起？

@jwezorek 感谢您的回复。它编译时没有任何警告。它在计算指数时挂起。它做了一些事情，但它从未计算过指数。即使将位数设置为 16（但仍使用 mpc_complex_backened）。

你能在调试器中单步执行吗？

是的。我试图尽可能地“进入”。调试器抛出错误“无法解析不存在的文件.../src/init2.c”。从 mpc_complex_imp() 调用函数 mpc_init2 时抛出此错误。可以跨过去。跨步时，它会挂起一个模板 BOOST_MP_CXX14_CONSTEXPR 编号（...）。也许有更好的方法来收集所有这些信息？我可以试试看。

有趣的是它在 init 函数中存在问题。你知道你是否安装了GMP？我不在 Linux 上，因此无法对此进行测试，但由于我没有 GMP，因此视觉工作室无法为我构建您的代码。也许您也没有，但它实际上可以与您的设置一起编译？

【参考方案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 位。