* Use the volatile qualifier for 'tiny' to prevent compilers

  from erronously constant folding expressions of the form
  '1 - tiny'.  This allows erf[f](x) to raise INEXACT.

* Use 0.5, 1, and 2, which are exactly representable in radix-2
  floating point formats.  This reduces diffs between s_erf[fl].c.

* While here, add a comment about efx and efx8.

2 files changed, 21 insertions, 10 deletions

diff --git a/lib/msun/src/s_erf.c b/lib/msun/src/s_erf.c
index 337f2e8..910951b 100644
--- a/lib/msun/src/s_erf.c
+++ b/lib/msun/src/s_erf.c
@@ -111,13 +111,20 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
+/* XXX Prevent compilers from erroneously constant folding: */
+static const volatile double tiny= 1e-300;
+
 static const double
-tiny	    = 1e-300,
-half=  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
-one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
-two =  2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
+half= 0.5,
+one = 1,
+two = 2,
 /* c = (float)0.84506291151 */
 erx =  8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
+/*
+ * In the domain [0, 2**-28], only the first term in the power series
+ * expansion of erf(x) is used.  The magnitude of the first neglected
+ * terms is less than 2**-84.
+ */
 efx =  1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
 efx8=  1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
 /*
diff --git a/lib/msun/src/s_erff.c b/lib/msun/src/s_erff.c
index b97ca1d..6dfdf87 100644
--- a/lib/msun/src/s_erff.c
+++ b/lib/msun/src/s_erff.c
@@ -19,13 +19,18 @@ __FBSDID("$FreeBSD$");
 #include "math.h"
 #include "math_private.h"
 
+/* XXX Prevent compilers from erroneously constant folding: */
+static const volatile float tiny = 1e-30;
+
 static const float
-tiny	    = 1e-30,
-half=  5.0000000000e-01, /* 0x3F000000 */
-one =  1.0000000000e+00, /* 0x3F800000 */
-two =  2.0000000000e+00, /* 0x40000000 */
+half= 0.5,
+one = 1,
+two = 2,
+erx = 8.42697144e-01,			/* 0x3f57bb00 */
 /*
- * Coefficients for approximation to erf on [0,0.84375]
+ * In the domain [0, 2**-14], only the first term in the power series
+ * expansion of erf(x) is used.  The magnitude of the first neglected
+ * terms is less than 2**-42.
  */
 efx =  1.2837916613e-01, /* 0x3e0375d4 */
 efx8=  1.0270333290e+00, /* 0x3f8375d4 */
@@ -43,7 +48,6 @@ qq3  = -1.98859419e-03F, /* -0x1.04a626p-9 */
  * Domain [0.84375, 1.25], range ~[-1.953e-11,1.940e-11]:
  * |(erf(x) - erx) - p(x)/q(x)| < 2**-36.
  */
-erx  =  8.42697144e-01F, /*  0x1.af7600p-1.  erf(1) rounded to 16 bits. */
 pa0  =  3.64939137e-06F, /*  0x1.e9d022p-19 */
 pa1  =  4.15109694e-01F, /*  0x1.a91284p-2 */
 pa2  = -1.65179938e-01F, /* -0x1.5249dcp-3 */