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
hvdijk committed May 27, 2024
1 parent d3e3bda commit 04fc3ba
Show file tree
Hide file tree
Showing 4 changed files with 51 additions and 45 deletions.
25 changes: 15 additions & 10 deletions test_conformance/math_brute_force/binary_double.cpp
View file
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
1 change: 1 addition & 0 deletions test_conformance/math_brute_force/function_list.h
View file
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
66 changes: 33 additions & 33 deletions test_conformance/math_brute_force/reference_math.cpp
View file
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_copysignl(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_copysignl(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_copysignl(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_copysignl(
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_copysignl(1.0, z); // really should just be an xor
z = reference_fabs(z); // remove the sign again

// reduce once more
Expand Down Expand Up @@ -1084,7 +1084,7 @@ float reference_copysign(float x, float y)
}


double reference_copysignd(double x, double y)
double reference_copysignl(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_copysignl(0.0, x);

if (absx < HEX_DBL(+, 1, 0, +, 53))
x = reference_trunc(x + reference_copysignd(0.5, x));
x = reference_trunc(x + reference_copysignl(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_copysignl((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_copysignl(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_copysignl(HEX_DBL(+, 1, 0, +, 52), x);
double rounded = (x + magic) - magic;
x = reference_copysignd(rounded, x);
x = reference_copysignl(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_copysignl(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_copysignl(0.5, x);
x = reference_fabs(x);
if (x > 1.0) return cl_make_nan();

Expand Down Expand Up @@ -2156,7 +2156,7 @@ static long double reduce1l(long double x)
}

// Find the nearest multiple of 2
const long double r = reference_copysignl(unit_exp, x);
const long double r = reference_copysignll(unit_exp, x);
long double z = x + r;
z -= r;

Expand Down Expand Up @@ -2983,7 +2983,7 @@ long double reference_rootnl(long double x, int i)

if (isinf(x))
{
if (i < 0) return reference_copysignl(0.0L, x);
if (i < 0) return reference_copysignll(0.0L, x);

return x;
}
Expand Down Expand Up @@ -3016,7 +3016,7 @@ long double reference_rootnl(long double x, int i)
DivideDD(&iHi, &iLo, 1.0, i);
x = reference_powl(x, iHi) * reference_powl(x, iLo);

return reference_copysignl(x, sign);
return reference_copysignll(x, sign);
}

long double reference_rsqrtl(long double x) { return 1.0L / sqrtl(x); }
Expand All @@ -3032,23 +3032,23 @@ long double reference_sinpil(long double x)
r = 1.0L - r;

// sinPi zeros have the same sign as x
if (r == 0.0L) return reference_copysignl(0.0L, x);
if (r == 0.0L) return reference_copysignll(0.0L, x);

return reference_sinl(r * M_PIL);
}

long double reference_tanpil(long double x)
{
// set aside the sign (allows us to preserve sign of -0)
long double sign = reference_copysignl(1.0L, x);
long double sign = reference_copysignll(1.0L, x);
long double z = reference_fabsl(x);

// if big and even -- caution: only works if x only has single precision
if (z >= HEX_LDBL(+, 1, 0, +, 53))
{
if (z == INFINITY) return x - x; // nan

return reference_copysignl(
return reference_copysignll(
0.0L, x); // tanpi ( n ) is copysign( 0.0, n) for even integers n.
}

Expand All @@ -3064,7 +3064,7 @@ long double reference_tanpil(long double x)
if ((i & 1) && z == 0.0L) sign = -sign;

// track changes to the sign
sign *= reference_copysignl(1.0L, z); // really should just be an xor
sign *= reference_copysignll(1.0L, z); // really should just be an xor
z = reference_fabsl(z); // remove the sign again

// reduce once more
Expand Down Expand Up @@ -3223,11 +3223,11 @@ int reference_signbitl(long double x) { return 0 != signbitl(x); }
#else
int reference_signbitl(long double x) { return 0 != signbit(x); }
#endif
long double reference_copysignl(long double x, long double y);
long double reference_copysignll(long double x, long double y);
long double reference_roundl(long double x);
long double reference_cbrtl(long double x);

long double reference_copysignl(long double x, long double y)
long double reference_copysignll(long double x, long double y)
{
// We hope that the long double to double conversion proceeds with sign
// fidelity, even for zeros and NaNs
Expand All @@ -3250,7 +3250,7 @@ long double reference_roundl(long double x)

#if defined(__MINGW32__) && defined(__x86_64__)
long double absx = reference_fabsl(x);
if (absx < 0.5L) return reference_copysignl(0.0L, x);
if (absx < 0.5L) return reference_copysignll(0.0L, x);
#endif
return round((double)x);
}
Expand Down Expand Up @@ -3305,7 +3305,7 @@ long double reference_cbrtl(long double x)
powxy = reference_scalblnl(powxy, m);
}

return reference_copysignl(powxy, x);
return reference_copysignll(powxy, x);
}

long double reference_rintl(long double x)
Expand Down Expand Up @@ -3405,7 +3405,7 @@ long double reference_asinhl(long double x)
long double absx = reference_fabsl(x);
if (absx < cutoff) return x;

long double sign = reference_copysignl(1.0L, x);
long double sign = reference_copysignll(1.0L, x);

if (absx <= 4.0 / 3.0)
{
Expand Down Expand Up @@ -3438,7 +3438,7 @@ long double reference_atanhl(long double x)
*/
if (isnan(x)) return x + x;

long double signed_half = reference_copysignl(0.5L, x);
long double signed_half = reference_copysignll(0.5L, x);
x = reference_fabsl(x);
if (x > 1.0L) return cl_make_nan();

Expand Down Expand Up @@ -4969,9 +4969,9 @@ static long double reference_scalblnl(long double x, long n)

if (reference_isinf(x)) return x;

if (x == 0.0L || n < -2200) return reference_copysignl(0.0L, x);
if (x == 0.0L || n < -2200) return reference_copysignll(0.0L, x);

if (n > 2200) return reference_copysignl(INFINITY, x);
if (n > 2200) return reference_copysignll(INFINITY, x);

if (n < 0)
{
Expand Down Expand Up @@ -5020,8 +5020,8 @@ static long double reference_scalblnl(long double x, long n)
e = (int)((u.l & 0x7ff0000000000000LL) >> 52) - 1022;
}
e += n;
if (e >= 2047) return reference_copysignl(INFINITY, x);
if (e < -51) return reference_copysignl(0.0, x);
if (e >= 2047) return reference_copysignll(INFINITY, x);
if (e < -51) return reference_copysignll(0.0, x);
if (e <= 0)
{
bias += (e - 1);
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_copysignl(result, x) : result;
}

double reference_sqrt(double x) { return sqrt(x); }
Expand Down Expand Up @@ -5391,7 +5391,7 @@ long double reference_acosl(long double x)
{
// we handle the x < 0 case as pi - acos(|x|)

long double sign = reference_copysignl(1.0L, x);
long double sign = reference_copysignll(1.0L, x);
long double fabsx = reference_fabsl(x);
head -= head * sign; // x > 0 ? 0 : pi.hi
tail -= tail * sign; // x > 0 ? 0 : pi.low
Expand Down
4 changes: 2 additions & 2 deletions test_conformance/math_brute_force/reference_math.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,6 @@ 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);
double reference_exp10(double);
double reference_exp2(double x);
double reference_expm1(double x);
Expand Down Expand Up @@ -205,7 +204,8 @@ long double reference_acoshl(long double x);
long double reference_asinhl(long double x);
long double reference_atanhl(long double x);
long double reference_cbrtl(long double x);
long double reference_copysignl(long double x, long double y);
double reference_copysignl(double x, double y);
long double reference_copysignll(long double x, long double y);
long double reference_exp10l(long double);
long double reference_exp2l(long double x);
long double reference_expm1l(long double x);
Expand Down

0 comments on commit 04fc3ba

Please sign in to comment.