erf function gives wrong results with pgcpp - PGI 10.4

[Jacques Desfosses <jacques.desfosses@…>]

Hi experts,

I compiled the following code with pgcpp / PGI 10.4 on linux.

/opt/pgi-10.4/linux86/10.4/bin/pgcpp main.cpp -o erf.exe

#include <math.h>
#include <boost/math/special_functions/erf.hpp>
int main()
{
   double val(1.0);
   printf("BOOST : %-20.15E\n", boost::math::erf(val));
   printf("MATH_H: %-20.15E\n", erf(val));
}

Output is as follows:

BOOST : 8.368481544380342E-01
MATH_H: 8.427007929497149E-01

The value computed by BOOST is incorrect. The same code will give the correct value when compiled with VS2012 with Windows 7-64 bits OS.

My linux box is running SUSE Linux Enterprise Desktop 10 SP2 (x86_64)

Can you please help troubleshoot the problem? Other basic statistic functions will also return wrong values such as cdf and quantile of standard distribution.

Thanks a lot,

Jacques Desfosses

[anonymous]

I'm going to be away for the next few days, and don't have access to the PGI compiler anyway, in the mean time can you try running with BOOST_MATH_INSTRUMENT defined and also try with GCC as the compiler?

Thanks, John Maddock

[Jacques Desfosses <jacques.desfosses@…>]

Here is the output when -DBOOST_MATH_INSTRUMENT is added to the compile line.

../boost/math/special_functions/erf.hpp:1049 result_type = long double ../boost/math/special_functions/erf.hpp:1050 value_type = long double ../boost/math/special_functions/erf.hpp:1051 precision_type = boost::math::policies::digits2<53> ../boost/math/special_functions/erf.hpp:1071 tag_type = mpl_::int_<53> ../boost/math/special_functions/erf.hpp:177 53-bit precision erf_imp called ../boost/math/special_functions/erf.hpp:1049 result_type = long double ../boost/math/special_functions/erf.hpp:1050 value_type = long double ../boost/math/special_functions/erf.hpp:1051 precision_type = boost::math::policies::digits2<53> ../boost/math/special_functions/erf.hpp:1071 tag_type = mpl_::int_<53> ../boost/math/special_functions/erf.hpp:177 53-bit precision erf_imp called ../boost/math/special_functions/erf.hpp:1049 result_type = long double ../boost/math/special_functions/erf.hpp:1050 value_type = long double ../boost/math/special_functions/erf.hpp:1051 precision_type = boost::math::policies::digits2<53> ../boost/math/special_functions/erf.hpp:1071 tag_type = mpl_::int_<53> ../boost/math/special_functions/erf.hpp:177 53-bit precision erf_imp called ../boost/math/special_functions/erf.hpp:1049 result_type = long double ../boost/math/special_functions/erf.hpp:1050 value_type = long double ../boost/math/special_functions/erf.hpp:1051 precision_type = boost::math::policies::digits2<53> ../boost/math/special_functions/erf.hpp:1071 tag_type = mpl_::int_<53> ../boost/math/special_functions/erf.hpp:177 53-bit precision erf_imp called ../boost/math/special_functions/erf.hpp:1049 result_type = long double ../boost/math/special_functions/erf.hpp:1050 value_type = long double ../boost/math/special_functions/erf.hpp:1051 precision_type = boost::math::policies::digits2<53> ../boost/math/special_functions/erf.hpp:1071 tag_type = mpl_::int_<53> ../boost/math/special_functions/erf.hpp:177 53-bit precision erf_imp called ../boost/math/special_functions/erf.hpp:1049 result_type = long double ../boost/math/special_functions/erf.hpp:1050 value_type = long double ../boost/math/special_functions/erf.hpp:1051 precision_type = boost::math::policies::digits2<53> ../boost/math/special_functions/erf.hpp:1071 tag_type = mpl_::int_<53> ../boost/math/special_functions/erf.hpp:177 53-bit precision erf_imp called ../boost/math/special_functions/erf.hpp:1049 result_type = double ../boost/math/special_functions/erf.hpp:1050 value_type = long double ../boost/math/special_functions/erf.hpp:1051 precision_type = boost::math::policies::digits2<53> ../boost/math/special_functions/erf.hpp:1071 tag_type = mpl_::int_<53> ../boost/math/special_functions/erf.hpp:177 53-bit precision erf_imp called BOOST : 8.506645778731253E-01 MATH_H: 8.427007929497149E-01

[Jacques Desfosses <jacques.desfosses@…>]

I also tried it out with g++ version 4.1.2 and obtained correct answers.

g++ (GCC) 4.1.2 20070115 (SUSE Linux)

BOOST : 8.427007929497149E-01 MATH_H: 8.427007929497149E-01

Unfortunately, there are compile errors with this compiler when using -DBOOST_MATH_INSTRUMENT so I am not able to provide the detailed output.

g++ -I../ main.cpp -DBOOST_MATH_INSTRUMENT -o erf_g++.exe In file included from ../boost/fusion/tuple/tuple.hpp:22,

from ../boost/fusion/tuple.hpp:10, from ../boost/fusion/include/tuple.hpp:10, from ../boost/math/tools/tuple.hpp:89, from ../boost/math/special_functions/detail/igamma_inverse.hpp:13, from ../boost/math/special_functions/gamma.hpp:1651, from ../boost/math/special_functions/erf.hpp:15, from main.cpp:2:

../boost/fusion/tuple/detail/preprocessed/tuple.hpp:21:7: warning: no newline at end of file ../boost/math/special_functions/fpclassify.hpp: In function â?~int boost::math::detail::fpclassify_imp(T, const boost::math::detail::generic_tag<true>&)â?T: ../boost/math/special_functions/fpclassify.hpp:132: error: â?~setprecisionâ?T is not a member of â?~stdâ?T ../boost/math/special_functions/fpclassify.hpp: In function â?~int boost::math::detail::fpclassify_imp(T, const boost::math::detail::generic_tag<false>&)â?T: ../boost/math/special_functions/fpclassify.hpp:176: error: â?~setprecisionâ?T is not a member of â?~stdâ?T ../boost/math/special_functions/fpclassify.hpp: In function â?~int boost::math::detail::fpclassify_imp(T, boost::math::detail::ieee_copy_all_bits_tag)â?T: ../boost/math/special_functions/fpclassify.hpp:186: error: â?~setprecisionâ?T is not a member of â?~stdâ?T ../boost/math/special_functions/fpclassify.hpp:190: error: â?~setprecisionâ?T is not a member of â?~stdâ?T ../boost/math/special_functions/fpclassify.hpp:192: error: â?~setprecisionâ?T is not a member of â?~stdâ?T ../boost/math/special_functions/fpclassify.hpp:193: error: â?~setprecisionâ?T is not a member of â?~stdâ?T ../boost/math/special_functions/fpclassify.hpp: In function â?~int boost::math::detail::fpclassify_imp(T, boost::math::detail::ieee_copy_leading_bits_tag)â?T: ../boost/math/special_functions/fpclassify.hpp:215: error: â?~setprecisionâ?T is not a member of â?~stdâ?T

[John Maddock]

Thanks for the extra information, I'm attaching an updated version of boost/math/special_functions/erf.hpp which prints some more debugging info. Can you give it a try and report back the output? The last few lines should be:

d:\data\boost\trunk\boost\math\special_functions\erf.hpp:177 53-bit precision erf_imp called
d:\data\boost\trunk\boost\math\special_functions\erf.hpp:269 Y = 0.40593576431274414
d:\data\boost\trunk\boost\math\special_functions\erf.hpp:270 P[0] = -0.098090592216281233
d:\data\boost\trunk\boost\math\special_functions\erf.hpp:271 Q[0] = 1
d:\data\boost\trunk\boost\math\special_functions\erf.hpp:273 result = 0.427583576155807
d:\data\boost\trunk\boost\math\special_functions\erf.hpp:275 result = 0.15729920705028513

Thanks, John.

[John Maddock]

Debugging version of erf.hpp

[Jacques Desfosses <jacques.desfosses@…>]

Hi John,

The output with the updated version of erf.hpp is below. I just noticed that when I compile with -DBOOST_MATH_INSTRUMENT the output from BOOST is still incorrect (8.506645778731253E-01) but different than the output from the version compiled without the define (8.368481544380342E-01).

Thanks!

Jacques

../boost/math/special_functions/erf.hpp:1055 result_type = long double ../boost/math/special_functions/erf.hpp:1056 value_type = long double ../boost/math/special_functions/erf.hpp:1057 precision_type = boost::math::policies::digits2<53> ../boost/math/special_functions/erf.hpp:1077 tag_type = mpl_::int_<53> ../boost/math/special_functions/erf.hpp:177 53-bit precision erf_imp called ../boost/math/special_functions/erf.hpp:1055 result_type = long double ../boost/math/special_functions/erf.hpp:1056 value_type = long double ../boost/math/special_functions/erf.hpp:1057 precision_type = boost::math::policies::digits2<53> ../boost/math/special_functions/erf.hpp:1077 tag_type = mpl_::int_<53> ../boost/math/special_functions/erf.hpp:177 53-bit precision erf_imp called ../boost/math/special_functions/erf.hpp:1055 result_type = long double ../boost/math/special_functions/erf.hpp:1056 value_type = long double ../boost/math/special_functions/erf.hpp:1057 precision_type = boost::math::policies::digits2<53> ../boost/math/special_functions/erf.hpp:1077 tag_type = mpl_::int_<53> ../boost/math/special_functions/erf.hpp:177 53-bit precision erf_imp called ../boost/math/special_functions/erf.hpp:269 Y = 0.40593576431274414 ../boost/math/special_functions/erf.hpp:270 P[0] = -0.098090592216281252 ../boost/math/special_functions/erf.hpp:271 Q[0] = -8.7793466483067508e+269 ../boost/math/special_functions/erf.hpp:272 z = 1.25 ../boost/math/special_functions/erf.hpp:274 result = 0.40593576431274414 ../boost/math/special_functions/erf.hpp:276 result = 0.068071006921468324 ../boost/math/special_functions/erf.hpp:1055 result_type = long double ../boost/math/special_functions/erf.hpp:1056 value_type = long double ../boost/math/special_functions/erf.hpp:1057 precision_type = boost::math::policies::digits2<53> ../boost/math/special_functions/erf.hpp:1077 tag_type = mpl_::int_<53> ../boost/math/special_functions/erf.hpp:177 53-bit precision erf_imp called ../boost/math/special_functions/erf.hpp:1055 result_type = long double ../boost/math/special_functions/erf.hpp:1056 value_type = long double ../boost/math/special_functions/erf.hpp:1057 precision_type = boost::math::policies::digits2<53> ../boost/math/special_functions/erf.hpp:1077 tag_type = mpl_::int_<53> ../boost/math/special_functions/erf.hpp:177 53-bit precision erf_imp called ../boost/math/special_functions/erf.hpp:1055 result_type = long double ../boost/math/special_functions/erf.hpp:1056 value_type = long double ../boost/math/special_functions/erf.hpp:1057 precision_type = boost::math::policies::digits2<53> ../boost/math/special_functions/erf.hpp:1077 tag_type = mpl_::int_<53> ../boost/math/special_functions/erf.hpp:177 53-bit precision erf_imp called ../boost/math/special_functions/erf.hpp:1055 result_type = double ../boost/math/special_functions/erf.hpp:1056 value_type = long double ../boost/math/special_functions/erf.hpp:1057 precision_type = boost::math::policies::digits2<53> ../boost/math/special_functions/erf.hpp:1077 tag_type = mpl_::int_<53> ../boost/math/special_functions/erf.hpp:177 53-bit precision erf_imp called ../boost/math/special_functions/erf.hpp:269 Y = 0.40593576431274414 ../boost/math/special_functions/erf.hpp:270 P[0] = -0.098090592216281252 ../boost/math/special_functions/erf.hpp:271 Q[0] = -8.7793466483067508e+269 ../boost/math/special_functions/erf.hpp:272 z = 1 ../boost/math/special_functions/erf.hpp:274 result = 0.40593576431274414 ../boost/math/special_functions/erf.hpp:276 result = 0.14933542212687465 BOOST : 8.506645778731253E-01 MATH_H: 8.427007929497149E-01

[John Maddock]

Thanks, the problem value is:

Q[0] = -8.7793466483067508e+269

With Q being initialized like this:

         static const T Q[] = {    
            1L,
            BOOST_MATH_BIG_CONSTANT(T, 53, 1.84759070983002217845),
            BOOST_MATH_BIG_CONSTANT(T, 53, 1.42628004845511324508),
            BOOST_MATH_BIG_CONSTANT(T, 53, 0.578052804889902404909),
            BOOST_MATH_BIG_CONSTANT(T, 53, 0.12385097467900864233),
            BOOST_MATH_BIG_CONSTANT(T, 53, 0.0113385233577001411017),
            BOOST_MATH_BIG_CONSTANT(T, 53, 0.337511472483094676155e-5),
         };

Can you take a look and see if the compiler has a problem with initialization from a long int? Initialization from long and long long is used all over the place in Boost.Math, so while I could easily patch that case, it would be hard to get all of them :-(

Besides it really should work!

[Jacques Desfosses <jacques.desfosses@…>]

Thanks again for your help! The compiler is okay with long initialization. The compiler seems to dislike the mixing of 1L and BOOST_MATH_BIG_CONSTANT initializers. I was able to obtain correct results by defining Q using these 2 methods:

1) Without BOOST_MATH_BIG_CONSTANT

static const T Q[] = {    
            1L,
            1.84759070983002217845,
            1.42628004845511324508,
            0.578052804889902404909,
            0.12385097467900864233,
            0.0113385233577001411017,
            0.337511472483094676155e-5,
         };

2) Only with BOOST_MATH_BIG_CONSTANT

static const T Q[] = {    
            BOOST_MATH_BIG_CONSTANT(T, 53, 1),
            BOOST_MATH_BIG_CONSTANT(T, 53, 1.84759070983002217845),
            BOOST_MATH_BIG_CONSTANT(T, 53, 1.42628004845511324508),
            BOOST_MATH_BIG_CONSTANT(T, 53, 0.578052804889902404909),
            BOOST_MATH_BIG_CONSTANT(T, 53, 0.12385097467900864233),
            BOOST_MATH_BIG_CONSTANT(T, 53, 0.0113385233577001411017),
            BOOST_MATH_BIG_CONSTANT(T, 53, 0.337511472483094676155e-5),
         };

I also had the opportunity to try out an older PGI compiler (version 7.1-3) and I also ran into the same problem. So the issue is not specific to PGI version 10.4.

[John Maddock]

(In [84714]) Don't mix literal and non-literal initializers in one table - it causes the PGI compiler to generate incorrect code. Refs #8621.

[John Maddock]

Can you try the patches referenced above?

Many thanks, John.

[Jacques Desfosses <jacques.desfosses@…>]

Hi John. Sorry for the delay. The patch works well. Did you, by any chance, forget to modify the following arrays?

boost/math/special_functions/detail/lgamma_small.hpp

static const T Q[] = {
   static_cast<T>(0.1e1),
   static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.196202987197795200688e1)),
   static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.148019669424231326694e1)),
   static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.541391432071720958364e0)),
   static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.988504251128010129477e-1)),
   static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.82130967464889339326e-2)),
   static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.224936291922115757597e-3)),
   static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.223352763208617092964e-6))
};


static const T Q[] = {
   static_cast<T>(0.1e1),
   static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.150169356054485044494e1)),
   static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.846973248876495016101e0)),
   static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.220095151814995745555e0)),
   static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.25582797155975869989e-1)),
   static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.100666795539143372762e-2)),
   static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.827193521891290553639e-6))
};

boost/math/special_functions/zeta.hpp

static const T Q[8] = {
   1,
   BOOST_MATH_BIG_CONSTANT(T, 64, 0.286577739726542730421),
   BOOST_MATH_BIG_CONSTANT(T, 64, 0.0447355811517733225843),
   BOOST_MATH_BIG_CONSTANT(T, 64, 0.00430125107610252363302),
   BOOST_MATH_BIG_CONSTANT(T, 64, 0.000284956969089786662045),
   BOOST_MATH_BIG_CONSTANT(T, 64, 0.116188101609848411329e-4),
   BOOST_MATH_BIG_CONSTANT(T, 64, 0.278090318191657278204e-6),
   BOOST_MATH_BIG_CONSTANT(T, 64, -0.19683620233222028478e-8),
};

I ran into another issue with the PGI compiler and the generic quantile finder. If it turns out to be an issue with boost, I will log another ticket. Thanks so much for your help,

Jacques

[John Maddock]

(In [84789]) Don't mix literal and non-literal initializers in one table - it causes the PGI compiler to generate incorrect code. Refs #8621.

[John Maddock]

(In [84864]) Apply patches from #8621. Fixes #8621.

[John Maddock]

(In [85987]) Merge accumulated patches from Trunk. Refs #8384, Refs #8855, refs #9107, refs #9109, refs #8333, refs #8621, refs #8732, refs #8733, refs #8837, refs #8940, refs #9042, refs #9087, refs #9104, refs #9126.