2

我是 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代码的 -type 版本,直到boost::multiprecision::mpfr_float_backend <2000>.

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

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

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

编译器:g++

cpp标准:gnu++17

什么可能是这个问题的根源?

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

4

1 回答 1

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 位。

于 2021-08-06T13:58:25.043 回答