Skip to content

Commit

Permalink
Disable implicit conversion of copysign arguments.
Browse files Browse the repository at this point in the history
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.
  • Loading branch information
hvdijk committed May 30, 2024
1 parent d3e3bda commit ca64b3e
Show file tree
Hide file tree
Showing 5 changed files with 49 additions and 29 deletions.
25 changes: 15 additions & 10 deletions test_conformance/math_brute_force/binary_double.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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();

Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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);

Expand Down Expand Up @@ -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 correct3 = func.f_ff(-0.0, s2[j]);
long double correct2 = ref_func(0.0, s2[j]);
long double correct3 = ref_func(-0.0, s2[j]);
float err2 =
Bruteforce_Ulp_Error_Double(test, correct2);
float err3 =
Expand All @@ -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);
correct3 = func.f_ff(-0.0, 0.0);
long double correct4 = func.f_ff(0.0, -0.0);
long double correct5 = func.f_ff(-0.0, -0.0);
correct2 = ref_func(0.0, 0.0);
correct3 = ref_func(-0.0, 0.0);
long double correct4 = ref_func(0.0, -0.0);
long double correct5 = ref_func(-0.0, -0.0);
err2 =
Bruteforce_Ulp_Error_Double(test, correct2);
err3 =
Expand Down Expand Up @@ -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 correct3 = func.f_ff(s[j], -0.0);
long double correct2 = ref_func(s[j], 0.0);
long double correct3 = ref_func(s[j], -0.0);
float err2 =
Bruteforce_Ulp_Error_Double(test, correct2);
float err3 =
Expand Down
16 changes: 15 additions & 1 deletion test_conformance/math_brute_force/function_list.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,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
Expand Down Expand Up @@ -212,7 +214,19 @@ const Func functionList[] = {
ENTRY(atan2pi, 6.0f, 6.0f, FTZ_OFF, binaryF),
ENTRY(cbrt, 2.0f, 4.0f, FTZ_OFF, unaryF),
ENTRY(ceil, 0.0f, 0.0f, FTZ_OFF, unaryF),
ENTRY(copysign, 0.0f, 0.0f, FTZ_OFF, binaryF),
{ "copysign",
"copysign",
{ (void*)reference_copysignf },
{ (void*)reference_copysign },
{ (void*)reference_copysignf },
0.0f,
0.0f,
0.0f,
INFINITY,
INFINITY,
FTZ_OFF,
RELAXED_OFF,
binaryF },
ENTRY_EXT(cos, 4.0f, 4.0f, 0.00048828125f, FTZ_OFF, unaryF,
0.00048828125f), // relaxed ulp 2^-11
ENTRY(cosh, 4.0f, 4.0f, FTZ_OFF, unaryF),
Expand Down
1 change: 1 addition & 0 deletions test_conformance/math_brute_force/function_list.h
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
32 changes: 16 additions & 16 deletions test_conformance/math_brute_force/reference_math.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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); }
Expand All @@ -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);
}
Expand All @@ -717,15 +717,15 @@ 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
if (z >= HEX_DBL(+, 1, 0, +, 24))
{
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.
}

Expand All @@ -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
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand All @@ -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;
}
Expand All @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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();

Expand Down Expand Up @@ -5290,7 +5290,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); }
Expand Down
4 changes: 2 additions & 2 deletions test_conformance/math_brute_force/reference_math.h
Original file line number Diff line number Diff line change
Expand Up @@ -87,8 +87,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);
double reference_copysignd(double x, double y);
float reference_copysignf(float x, float y);
double reference_copysign(double x, double y);
double reference_exp10(double);
double reference_exp2(double x);
double reference_expm1(double x);
Expand Down

0 comments on commit ca64b3e

Please sign in to comment.