Disable implicit conversion of copysign arguments. (#1970)

In binary_float.cpp, copysign is special cased. All the reference
functions there take double arguments, except for reference_copysign
which takes float. This commit copies that approach to special case
copysign in binary_double.cpp as well: all the reference functions there
take long double arguments, but this commit changes reference_copysignl
to take double. The rationale for this in binary_float applies equally
to binary_double: conversions of NAN are not required to preserve its
sign bit. On architectures where conversion of NAN resets the sign bit,
copysign fp64 would return incorrect reference results.

test_conformance/math_brute_force/binary_double.cpp

@@ -219,6 +219,7 @@ cl_int Test(cl_uint job_id, cl_uint thread_id, void *data)
     cl_double *r;
     cl_double *s;
     cl_double *s2;
+    cl_int copysign_test = 0;
 
     Force64BitFPUPrecision();
 
@@ -377,12 +378,16 @@ cl_int Test(cl_uint job_id, cl_uint thread_id, void *data)
 
     if (gSkipCorrectnessTesting) return CL_SUCCESS;
 
+    if (!strcmp(name, "copysign")) copysign_test = 1;
+
+#define ref_func(s, s2) (copysign_test ? func.f_ff_d(s, s2) : func.f_ff(s, s2))
+
     // Calculate the correctly rounded reference result
     r = (cl_double *)gOut_Ref + thread_id * buffer_elements;
     s = (cl_double *)gIn + thread_id * buffer_elements;
     s2 = (cl_double *)gIn2 + thread_id * buffer_elements;
     for (size_t j = 0; j < buffer_elements; j++)
-        r[j] = (cl_double)func.f_ff(s[j], s2[j]);
+        r[j] = (cl_double)ref_func(s[j], s2[j]);
 
     // Read the data back -- no need to wait for the first N-1 buffers but wait
     // for the last buffer. This is an in order queue.
@@ -412,7 +417,7 @@ cl_int Test(cl_uint job_id, cl_uint thread_id, void *data)
             if (t[j] != q[j])
             {
                 cl_double test = ((cl_double *)q)[j];
-                long double correct = func.f_ff(s[j], s2[j]);
+                long double correct = ref_func(s[j], s2[j]);
                 float err = Bruteforce_Ulp_Error_Double(test, correct);
                 int fail = !(fabsf(err) <= ulps);
 
@@ -449,8 +454,8 @@ cl_int Test(cl_uint job_id, cl_uint thread_id, void *data)
                         // retry per section 6.5.3.3
                         if (IsDoubleSubnormal(s[j]))
                         {
-                            long double correct2 = func.f_ff(0.0, s2[j]);
+                            long double correct2 = ref_func(0.0, s2[j]);
-                            long double correct3 = func.f_ff(-0.0, s2[j]);
+                            long double correct3 = ref_func(-0.0, s2[j]);
                             float err2 =
                                 Bruteforce_Ulp_Error_Double(test, correct2);
                             float err3 =
@@ -472,10 +477,10 @@ cl_int Test(cl_uint job_id, cl_uint thread_id, void *data)
                             // try with both args as zero
                             if (IsDoubleSubnormal(s2[j]))
                             {
-                                correct2 = func.f_ff(0.0, 0.0);
+                                correct2 = ref_func(0.0, 0.0);
-                                correct3 = func.f_ff(-0.0, 0.0);
+                                correct3 = ref_func(-0.0, 0.0);
-                                long double correct4 = func.f_ff(0.0, -0.0);
+                                long double correct4 = ref_func(0.0, -0.0);
-                                long double correct5 = func.f_ff(-0.0, -0.0);
+                                long double correct5 = ref_func(-0.0, -0.0);
                                 err2 =
                                     Bruteforce_Ulp_Error_Double(test, correct2);
                                 err3 =
@@ -507,8 +512,8 @@ cl_int Test(cl_uint job_id, cl_uint thread_id, void *data)
                         }
                         else if (IsDoubleSubnormal(s2[j]))
                         {
-                            long double correct2 = func.f_ff(s[j], 0.0);
+                            long double correct2 = ref_func(s[j], 0.0);
-                            long double correct3 = func.f_ff(s[j], -0.0);
+                            long double correct3 = ref_func(s[j], -0.0);
                             float err2 =
                                 Bruteforce_Ulp_Error_Double(test, correct2);
                             float err3 =

test_conformance/math_brute_force/function_list.cpp

@@ -73,6 +73,8 @@
 #define binaryF_two_results_i NULL
 #define mad_function NULL
 
+#define reference_copysignf NULL
+#define reference_copysign NULL
 #define reference_sqrt NULL
 #define reference_sqrtl NULL
 #define reference_divide NULL
@@ -250,7 +252,20 @@ const Func functionList[] = {
     ENTRY(atan2pi, 6.0f, 6.0f, 2.0f, FTZ_OFF, binaryF),
     ENTRY(cbrt, 2.0f, 4.0f, 2.f, FTZ_OFF, unaryF),
     ENTRY(ceil, 0.0f, 0.0f, 0.f, FTZ_OFF, unaryF),
-    ENTRY(copysign, 0.0f, 0.0f, 0.f, FTZ_OFF, binaryF),
+    { "copysign",
+      "copysign",
+      { (void*)reference_copysignf },
+      { (void*)reference_copysign },
+      { (void*)reference_copysignf },
+      0.0f,
+      0.0f,
+      0.0f,
+      0.0f,
+      INFINITY,
+      INFINITY,
+      FTZ_OFF,
+      RELAXED_OFF,
+      binaryF },
     ENTRY_EXT(cos, 4.0f, 4.0f, 2.f, 0.00048828125f, FTZ_OFF, unaryF,
               0.00048828125f), // relaxed ulp 2^-11
     ENTRY(cosh, 4.0f, 4.0f, 2.f, FTZ_OFF, unaryF),

test_conformance/math_brute_force/function_list.h

@@ -52,6 +52,7 @@ union dptr {
     long double (*f_f)(long double);
     long double (*f_u)(cl_ulong);
     int (*i_f)(long double);
+    double (*f_ff_d)(double, double);
     long double (*f_ff)(long double, long double);
     int (*i_ff)(long double, long double);
     long double (*f_fi)(long double, int);

test_conformance/math_brute_force/reference_math.cpp

@@ -691,7 +691,7 @@ double reference_rootn(double x, int i)
     double sign = x;
     x = reference_fabs(x);
     x = reference_exp2(reference_log2(x) / (double)i);
-    return reference_copysignd(x, sign);
+    return reference_copysign(x, sign);
 }
 
 double reference_rsqrt(double x) { return 1.0 / reference_sqrt(x); }
@@ -707,7 +707,7 @@ double reference_sinpi(double x)
         r = 1 - r;
 
     // sinPi zeros have the same sign as x
-    if (r == 0.0) return reference_copysignd(0.0, x);
+    if (r == 0.0) return reference_copysign(0.0, x);
 
     return reference_sin(r * M_PI);
 }
@@ -717,7 +717,7 @@ double reference_relaxed_sinpi(double x) { return reference_sinpi(x); }
 double reference_tanpi(double x)
 {
     // set aside the sign  (allows us to preserve sign of -0)
-    double sign = reference_copysignd(1.0, x);
+    double sign = reference_copysign(1.0, x);
     double z = reference_fabs(x);
 
     // if big and even  -- caution: only works if x only has single precision
@@ -725,7 +725,7 @@ double reference_tanpi(double x)
     {
         if (z == INFINITY) return x - x; // nan
 
-        return reference_copysignd(
+        return reference_copysign(
             0.0, x); // tanpi ( n ) is copysign( 0.0, n)  for even integers n.
     }
 
@@ -739,7 +739,7 @@ double reference_tanpi(double x)
     if ((i & 1) && z == 0.0) sign = -sign;
 
     // track changes to the sign
-    sign *= reference_copysignd(1.0, z); // really should just be an xor
+    sign *= reference_copysign(1.0, z); // really should just be an xor
     z = reference_fabs(z); // remove the sign again
 
     // reduce once more
@@ -1070,7 +1070,7 @@ int reference_signbit(float x) { return 0 != signbit(x); }
 // Missing functions for win32
 
 
-float reference_copysign(float x, float y)
+float reference_copysignf(float x, float y)
 {
     union {
         float f;
@@ -1084,7 +1084,7 @@ float reference_copysign(float x, float y)
 }
 
 
-double reference_copysignd(double x, double y)
+double reference_copysign(double x, double y)
 {
     union {
         double f;
@@ -1101,10 +1101,10 @@ double reference_copysignd(double x, double y)
 double reference_round(double x)
 {
     double absx = reference_fabs(x);
-    if (absx < 0.5) return reference_copysignd(0.0, x);
+    if (absx < 0.5) return reference_copysign(0.0, x);
 
     if (absx < HEX_DBL(+, 1, 0, +, 53))
-        x = reference_trunc(x + reference_copysignd(0.5, x));
+        x = reference_trunc(x + reference_copysign(0.5, x));
 
     return x;
 }
@@ -1115,7 +1115,7 @@ double reference_trunc(double x)
     {
         cl_long l = (cl_long)x;
 
-        return reference_copysignd((double)l, x);
+        return reference_copysign((double)l, x);
     }
 
     return x;
@@ -1132,16 +1132,16 @@ double reference_trunc(double x)
 
 double reference_cbrt(double x)
 {
-    return reference_copysignd(reference_pow(reference_fabs(x), 1.0 / 3.0), x);
+    return reference_copysign(reference_pow(reference_fabs(x), 1.0 / 3.0), x);
 }
 
 double reference_rint(double x)
 {
     if (reference_fabs(x) < HEX_DBL(+, 1, 0, +, 52))
     {
-        double magic = reference_copysignd(HEX_DBL(+, 1, 0, +, 52), x);
+        double magic = reference_copysign(HEX_DBL(+, 1, 0, +, 52), x);
         double rounded = (x + magic) - magic;
-        x = reference_copysignd(rounded, x);
+        x = reference_copysign(rounded, x);
     }
 
     return x;
@@ -1174,7 +1174,7 @@ double reference_asinh(double x)
     double absx = reference_fabs(x);
     if (absx < HEX_DBL(+, 1, 0, -, 28)) return x;
 
-    double sign = reference_copysignd(1.0, x);
+    double sign = reference_copysign(1.0, x);
 
     if (absx > HEX_DBL(+, 1, 0, +, 28))
         return sign
@@ -1206,7 +1206,7 @@ double reference_atanh(double x)
      */
     if (isnan(x)) return x + x;
 
-    double signed_half = reference_copysignd(0.5, x);
+    double signed_half = reference_copysign(0.5, x);
     x = reference_fabs(x);
     if (x > 1.0) return cl_make_nan();
 
@@ -5333,7 +5333,7 @@ double reference_pow(double x, double y)
     __log2_ep(&hi, &lo, fabsx);
     double prod = y * hi;
     double result = reference_exp2(prod);
-    return isOddInt ? reference_copysignd(result, x) : result;
+    return isOddInt ? reference_copysign(result, x) : result;
 }
 
 double reference_sqrt(double x) { return sqrt(x); }

test_conformance/math_brute_force/reference_math.h

@@ -88,8 +88,8 @@ double reference_acosh(double x);
 double reference_asinh(double x);
 double reference_atanh(double x);
 double reference_cbrt(double x);
-float reference_copysign(float x, float y);
+float reference_copysignf(float x, float y);
-double reference_copysignd(double x, double y);
+double reference_copysign(double x, double y);
 double reference_exp10(double);
 double reference_exp2(double x);
 double reference_expm1(double x);