Modify

Opened 4 years ago

Closed 4 years ago

Last modified 4 years ago

#8621 closed Bugs (fixed)

erf function gives wrong results with pgcpp - PGI 10.4

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

Description (last modified by viboes)

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 (1)

erf.hpp (53.0 KB) - added by John Maddock 4 years ago.
Debugging version of erf.hpp

Download all attachments as: .zip

Change History (16)

comment:1 Changed 4 years ago by viboes

Description: modified (diff)

comment:2 Changed 4 years ago by viboes

Component: Nonemath
Owner: set to John Maddock

comment:3 Changed 4 years 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 4 years 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 4 years 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 4 years ago by 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.

Changed 4 years ago by John Maddock

Attachment: erf.hpp added

Debugging version of erf.hpp

comment:7 Changed 4 years 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 4 years ago by 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!

comment:9 Changed 4 years 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 4 years ago by 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.

comment:11 Changed 4 years ago by John Maddock

Can you try the patches referenced above?

Many thanks, John.

comment:12 Changed 4 years 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 4 years ago by 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.

comment:14 Changed 4 years ago by John Maddock

Resolution: fixed
Status: newclosed

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

comment:15 Changed 4 years ago by 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.

Modify Ticket

Change Properties
Set your email in Preferences
Action
as closed The owner will remain John Maddock.
The resolution will be deleted.

Add Comment


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

 
Note: See TracTickets for help on using tickets.