Re: [boost][multiprecision/mpc][math/quadrature] Why boost::multiprecision::exp gets stuck when evaluating complex-valued integral?

From: Kirill Kudashkin (kirill.kudashkin_at_[hidden])
Date: 2021-08-06 10:18:08

thank you for your reply. I implemented your idea and it worked.

Though, I had to constraint the exponent not only from below but from above.

Thanks and

Cheers,

Kirill

On 8/5/21 6:54 PM, John Maddock via Boost-users wrote:
>
>> I have encountered a problem while using Boost C++ libraries ( 
>> multiprecision/mpc.hpp, math/quadrature) which I cannot resolve.
>>
>> This problems appears when running the following piece of code
>>
>> ||
>
> OK, the value that triggers the issue is 
> -6.275135555002646369320950230099689165841342351751646e-118909138 - 
> 6.275135555002646369320950230099689165841342351751646e-118909138i
>
> And this appears to be an exponent range issue - basically mpc_exp 
> gets slower and slower as the exponents approach -INF and eventually 
> grinds to a halt.
>
> I'm not yet sure if this is a bug in MPC, or if you've just hit a 
> stupidly hard value to compute - wolframalpha dies as well with this 
> case.
>
> For the purposes of calculating the integral, guess you could just 
> return 1 from your functor whenever the exponent is sufficiently small?
>
> HTH, John.
>
>> |#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; } |
>>
>> ||Namely, the problem appears when evaluating the exponent.
>>
>> The code is compiled without any problems, but hangs on computing
>>
>> |Â Â Â  mpc_type ex2 = mpns::exp(-z*x2) ;|
>>
>> It was possible to go into the function while debugging.
>>
>> It crashed eventually, when the compiler reached "unable to resolve 
>> non-existing file.../src/init2.c".
>>
>> I let it run for a while checking the call stack. It appears that the 
>> programs is stuck atcalling ___gmpn_sub_n_coreisbr().
>>
>> Or better, in the call stuck it appeared frequently.
>>
>>
>> I did several tests.
>>
>> The real-valued version (with mpfr_backend <2000>) works without any 
>> problems.
>>
>> cpp_complex_backend <100> works fine, but I could not test higher 
>> precision yet, since it is much slower.
>>
>> I found a work-around which is not ideal. Namely,
>>
>> mpc_type ex2 = mpc_type(mpns::exp(x2)) ;
>> Â Â Â  mpc_type ezx2 = mpns::pow(ex2,-z) ;
>>
>> works fine with mpc_complex_backend <2000>.
>>
>> The lambda-function as it is declared in the code above works and 
>> yields correct numbers.
>>
>> The problem appears when it is called from the integrator.
>>
>>
>> My system.
>>
>> I use openSUSE Tumbeleweed.
>>
>> I installed boost libraries as provided by 
>> |boost_1_76_0-gnu-mpich-hpc-devel.|
>>
>> Complier -> g++, cppStandard gnu+17.
>>
>> I check that I have all libs, e.g. libgmp, libmpc, libmpfr and etc.
>>
>> All headers are present (the programs is complied without any warnings).
>>
>>
>> Thank you in advance and
>>
>> Best wishes,
>>
>> Kirill
>>
>> ||
>>
>>
>>
>> _______________________________________________
>> Boost-users mailing list
>> Boost-users_at_[hidden]
>> https://lists.boost.org/mailman/listinfo.cgi/boost-users
>
>
>