STRICT_ASSIGN macro is broken

[kargl]

Testing openlibm's single argument float functions on i686-*-freebsd has revealed a significant problem with accuracy; especially for exp2f(). In the following, ULP1 is the max ULP recorded from exhaustive testing in the indicated interval when openlibm is built without change on i686-*-freebsd. ULP2 is recorded when STRICT_ASSIGN() in math_private.h is restored its original form obtained from FreeBSD. It is unclear why STRICT_ASSIGN() was changed. The starred (*) items in the table show the importance of volatile in STRICT_ASSIGN().

function   interval        ULP1      ULP2
acosf       [0.5,1]        0.61275   0.61275
acoshf      [1,2]          2.61417   1.60654  *
asinf       [0.5,1]        0.69603   0.69603
asinhf      [1,2]          1.87656   1.08336  *
atanf       [0.5,10]       0.51169   0.51169
atanhf      [0.5,0.95]     1.24099   0.89150  *
cbrtf       [1,10]         0.50000   0.50000
cosf        [1,3.14159]    0.50089   0.50089
coshf       [1,10]         0.75235   0.75248
erff        [0.5,2]        0.55782   0.55782
erfcf       [0.5,2]        1.54599   1.54599
expf        [1,2]          0.50915   0.53298  *
exp2f       [1,8]          359512.0  22729.4  *
expm1f      [0.0625,1]     0.50013   0.50668
j0f         [1,10]         2945366.  2945366.
j1f         [1,10]         9107852.  9107852.
lgammaf     [1,3]          0.98491   0.98491
logf        [1,2]          0.50000   0.50000
log2f       [1,8]          0.50042   0.50042
log1pf      [0.5,1.5]      1.83431   0.50196
log10f      [1,8]          0.50000   0.50000
sinf        [1,3.14159]    0.50089   0.50089
sinhf       [1,10]         0.75246   0.75207  *
sqrtf       [1,10]         0.50000   0.50000
tanf        [0.0625,0.785] 0.79979   0.79979
tanhf       [1,10]         0.51103   0.51103
tgammaf     [1,3]          0.50000   0.50000
y0f         [1,10]         4491818.  4491818.
y1f         [1,10]         2956156.  2956156.

A deeper inspection of exp2f() shows the following results for the original unmodified openlibm.

% ./tlibm_olibm -PDEf exp2
Interval tested for exp2f: [1,8]
       ulp <= 0.5:  0.936%    235635 |   0.936%    235635
0.5 <  ulp <  0.6:  0.016%      4098 |   0.953%    239733
0.6 <  ulp <  0.7:  0.024%      6147 |   0.977%    245880
0.7 <  ulp <  0.8:  0.008%      2049 |   0.985%    247929
0.8 <  ulp <  0.9:  0.016%      4098 |   1.001%    252027
0.9 <  ulp <  1.0:  0.008%      2049 |   1.010%    254076
1.0 <  ulp <  1.5:  0.334%     84009 |   1.343%    338085
1.5 <  ulp <  2.0:  0.187%     47127 |   1.531%    385212
2.0 <  ulp <  3.0:  0.839%    211047 |   2.369%    596259
3.0 <  ulp       : 97.631%  24569565 | 100.000%  25165824
Max ulp: 359512.000000 at 1.96875286e+00

After restoring STRICT_ASSIGN(), one has

% ./tlibm_olibm -PDEf exp2
Interval tested for exp2f: [1,8]
       ulp <= 0.5:  0.056%     14072 |   0.056%     14072
0.5 <  ulp <  0.6:  0.000%         8 |   0.056%     14080
3.0 <  ulp <     : 99.944%  25151744 | 100.000%  25165824
Max ulp: 22729.386719 at 1.00195301e+00

Strictly (pun intended) speaking, this is an improvement. But, comparison exp2f() from FreeBSD libm shows

% ./tlibm_libm -PDEf exp2
Interval tested for exp2f: [1,8]
       ulp <= 0.5: 99.959%  25155610 |  99.959%  25155610
0.5 <  ulp <  0.6:  0.041%     10214 | 100.000%  25165824
Max ulp: 0.500980 at 1.97115958e+00

So, clearly openlibm is somewhat broken on i686-*-freebsd. This is attributed to the broken STRICT_ASSIGN(), and the use of -march=i686. If I change Make.inc to use -march=core2 or my CPU, then openlibm agrees with FreeBSD libm. A better method of choosing the -march seems to be desired.

[kargl]

Here's the patch that was used to restore STRICT_ASSIGN().

diff --git a/src/math_private.h b/src/math_private.h
index 17027d6..dcfe96f 100644
--- a/src/math_private.h
+++ b/src/math_private.h
@@ -203,10 +203,10 @@ do {                              \
 } while (0)

+#ifndef __FreeBSD__
 //VBS
 #define    STRICT_ASSIGN(type, lval, rval) ((lval) = (rval))
-
-/* VBS
+#else
 #ifdef FLT_EVAL_METHOD
 // Attempt to get strict C99 semantics for assignment with non-C99 compilers.
 #if FLT_EVAL_METHOD == 0 || __GNUC__ == 0
@@ -215,7 +215,7 @@ do {                                \
 #define    STRICT_ASSIGN(type, lval, rval) do {    \
    volatile type __lval;           \
                        \
-   if (sizeof(type) >= sizeof(double)) \
+   if (sizeof(type) >= sizeof(long double))    \
        (lval) = (rval);        \
    else {                  \
        __lval = (rval);        \
@@ -224,7 +224,7 @@ do {                                \
 } while (0)
 #endif
 #endif
-*/
+#endif

 /*
  * Common routine to process the arguments to nan(), nanf(), and nanl().

[kargl]

Having now tested exp2(x), it is clear that the STRICT_ASSIGN() macro in math_private.h should be restored to its definition when inherited from FreeBSD. It is needed for proper rounding. Testing 5 million values in the interval [1,8] with FreeBSD libm returns

./tlibm_libm -d -N 5 exp2
Interval tested for exp2: [1,8]
count: 5000000
  xm =  5.7284883456976692e+00, /* 0x4016e9f8, 0xd951278f */
libm =  5.3020866609967030e+01, /* 0x404a82ab, 0xc1cfb182 */
mpfr =  5.3020866609967037e+01, /* 0x404a82ab, 0xc1cfb183 */
 ULP = 0.50264

The result with OpenLibm with no changes gives

 ./tlibm_olibm -d -N 5 exp2
Interval tested for exp2: [1,8]
count: 5000000
  xm =  6.9980477996095596e+00, /* 0x401bfe00, 0x3e0cabd9 */
libm =  1.2800008203401148e+02, /* 0x40600000, 0xac09acb0 */
mpfr =  1.2782691237307291e+02, /* 0x405ff4ec, 0x21dfc061 */
 ULP = 12185731569919.41602

[ViralBShah]

@kargl Can you open a PR here?

[ViralBShah]

It would be great to get the accumulated improvements from msun back into openlibm, but it is far from easy to do something like that.

[kargl]

Viral, I'm not sure what you mean by opening a PR? Isn't this issues/215 a PR? There is a patch in my Sep 8, 2020 comment. [edited] PR means Problem Report to me. Do you mean Pull Request? Have no idea how to do that.

[ViralBShah]

Haha - yes, I meant a Pull Request. I can convert your patch into a PR if you are not familiar with the process on github.