Unable to produce proper multiprecision output from sph_neumann() function

[James_Booth@…]

I am trying to use the boost c++ libraries to compute a multiprecision value for the spherical bessel functions and spherical neumann functions. To test the boost package I am using the relationships:

pi/6*Spherical_Bessel(0,pi/6) = sin(pi/6) = 0.5

-pi/6*Spherical_Neumann(0,pi/6) = cos(pi/6) = sqrt(3)/2 = 0.866.... My test code is:

#include <boost/multiprecision/mpfr.hpp>

#include <boost/math/special_functions/bessel.hpp>

#include <boost/math/constants/constants.hpp>

#include <boost/multiprecision/number.hpp>

#define PRECISION 40 Can be set to any number desired.

using namespace std; using namespace boost::multiprecision; using namespace boost;

mpfr_float::default_precision(PRECISION);

mpfr_float z, my_pi;

my_pi = boost::multiprecision::atan(mpfr_float("1.0"))*mpfr_float("4.0"); Define pi as a multiprecision value to precision = PRECISION.

z = my_pi/mpfr_float("6.0");

cout << z*mpfr_float(boost::math::sph_bessel(0,z)) << endl;

cout << mpfr_float("-1.0")*z*mpfr_float(boost::math::sph_neumann(0,z)) << endl;

========================================================================

The value of my_pi agrees with the value computed by Wolfram (to at least 100 decimal places). The value of z also agrees with the value from Wolfram. The output for the z*sph_bessel(0,z) function is correct: 0.5000000000.... The output for the -z*sph_neumann(0,z) function is incorrect: 0.86602540378443864*30*... instead of 0.86602540378443864*67*.... Which looks like the answer is not multiprecision but a float.

Can you please help with this bug? Thanks.

[John Maddock]

Confirmed - it's an issue with using mpfr_float (a variable precision type and no numeric_limits support) vs a fixed precision type like mpfr_float_50 which does produce the correct result. There's at least one bug somewhere, but I can't find it at present.

[John Maddock]

This required rather a lot of changes to fix - but is now fixed in develop - you will need both Math and Multiprecision lib develop branches though to see the fix. On the plus side you can now change precision in mid program and have it all still work... something we didn't support before.