Modify

Ticket #8621 (closed Bugs: fixed)

Opened 11 months ago

Last modified 7 months ago

erf function gives wrong results with pgcpp - PGI 10.4

Reported by: Jacques Desfosses <jacques.desfosses@…> Owned by: johnmaddock
Milestone: To Be Determined Component: math
Version: Boost 1.53.0 Severity: Problem
Keywords: erf pgi Cc:

Description (last modified by viboes) (diff)

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

Attachments

erf.hpp Download (53.0 KB) - added by johnmaddock 10 months ago.
Debugging version of erf.hpp

Change History

comment:1 Changed 11 months ago by viboes

  • Description modified (diff)

comment:2 Changed 11 months ago by viboes

  • Owner set to johnmaddock
  • Component changed from None to math

comment:3 Changed 11 months ago by 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

comment:4 Changed 11 months ago by 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

comment:5 Changed 11 months ago by 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

comment:6 Changed 10 months ago by johnmaddock

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.

Changed 10 months ago by johnmaddock

Debugging version of erf.hpp

comment:7 Changed 10 months ago by 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

comment:8 Changed 10 months ago by johnmaddock

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!

comment:9 Changed 10 months ago by 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.

comment:10 Changed 10 months ago by johnmaddock

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

comment:11 Changed 10 months ago by johnmaddock

Can you try the patches referenced above?

Many thanks, John.

comment:12 Changed 10 months ago by 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

comment:13 Changed 10 months ago by johnmaddock

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

comment:14 Changed 10 months ago by johnmaddock

  • Status changed from new to closed
  • Resolution set to fixed

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

comment:15 Changed 7 months ago by johnmaddock

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

View

Add a comment

Modify Ticket

Change Properties
<Author field>
Action
as closed
The resolution will be deleted. Next status will be 'reopened'
Author


E-mail address and user name can be saved in the Preferences.

 
Note: See TracTickets for help on using tickets.