Boost C++ Libraries: Ticket #8621: erf function gives wrong results with pgcpp - PGI 10.4

description changed

viboes — Mon, 03 Jun 2013 18:03:00 GMT

description modified (diff)

component changed; owner set

viboes — Mon, 03 Jun 2013 18:03:22 GMT

owner set to John Maddock
component None → math

anonymous — Mon, 03 Jun 2013 18:24:16 GMT

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

Mon, 03 Jun 2013 18:59:59 GMT

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

Mon, 03 Jun 2013 19:09:55 GMT

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 — Fri, 07 Jun 2013 18:08:04 GMT

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.

attachment set

John Maddock — Fri, 07 Jun 2013 18:11:59 GMT

attachment → erf.hpp

Debugging version of erf.hpp

Fri, 07 Jun 2013 18:26:16 GMT

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 — Sat, 08 Jun 2013 11:04:53 GMT

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!

Sun, 09 Jun 2013 03:37:52 GMT

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 — Sun, 09 Jun 2013 14:14:12 GMT

(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 — Sun, 09 Jun 2013 14:14:57 GMT

Can you try the patches referenced above?

Many thanks, John.

Thu, 13 Jun 2013 17:41:47 GMT

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 — Sat, 15 Jun 2013 08:04:56 GMT

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

status changed; resolution set

John Maddock — Fri, 21 Jun 2013 16:08:11 GMT

status new → closed
resolution → fixed

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

John Maddock — Sat, 28 Sep 2013 16:19:56 GMT

(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.