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

Dear all, 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 || |#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 ||

I have encountered a problem while using Boost C++ libraries ( multiprecision/mpc.hpp, math/quadrature) which I cannot resolve.
I'm a bit tied up for a few days, but I'll look into this as soon as I can, John.
This problems appears when running the following piece of code
|| |#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@lists.boost.org https://lists.boost.org/mailman/listinfo.cgi/boost-users
-- This email has been checked for viruses by Avast antivirus software. https://www.avast.com/antivirus

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@lists.boost.org https://lists.boost.org/mailman/listinfo.cgi/boost-users
-- This email has been checked for viruses by Avast antivirus software. https://www.avast.com/antivirus

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@lists.boost.org https://lists.boost.org/mailman/listinfo.cgi/boost-users

On 06/08/2021 11:18, Kirill Kudashkin via Boost-users wrote:
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.
Nod. One of the problems with this form of quadrature is that it can reach very far out into the end points.... sometimes you really need this though if that's where all the area is. -- This email has been checked for viruses by Avast antivirus software. https://www.avast.com/antivirus
participants (2)
-
John Maddock
-
Kirill Kudashkin