Enabled tgamma test for bruteforce suite (#2251)
Fixes #507 according to work plan from issue description.
This commit is contained in:
Marcin Hajder
GitHub
parent
4a923c074b
commit
3f4bad4a77
7 changed files with 323 additions and 15 deletions
|
|
@ -405,7 +405,7 @@ const Func functionList[] = {
|
|||
|
||||
ENTRY(tanh, 5.0f, 5.0f, 2.0f, 3.0f, FTZ_OFF, unaryF),
|
||||
ENTRY(tanpi, 6.0f, 6.0f, 2.0f, 3.0f, FTZ_OFF, unaryF),
|
||||
//ENTRY(tgamma, 16.0f, 16.0f, FTZ_OFF, unaryF), Commented this out until we can be sure this requirement is realistic
|
||||
ENTRY(tgamma, 16.0f, 16.0f, 4.f, 4.f, FTZ_OFF, unaryF),
|
||||
ENTRY(trunc, 0.0f, 0.0f, 0.0f, 0.0f, FTZ_OFF, unaryF),
|
||||
|
||||
HALF_ENTRY(cos, 8192.0f, 8192.0f, FTZ_ON, unaryOF),
|
||||
|
|
|
|||
|
|
@ -342,6 +342,7 @@ DO_TEST(sqrt_cr)
|
|||
DO_TEST(tan)
|
||||
DO_TEST(tanh)
|
||||
DO_TEST(tanpi)
|
||||
DO_TEST(tgamma)
|
||||
DO_TEST(trunc)
|
||||
DO_TEST(half_cos)
|
||||
DO_TEST(half_divide)
|
||||
|
|
|
|||
|
|
@ -2119,6 +2119,268 @@ int reference_not(double x)
|
|||
return r;
|
||||
}
|
||||
|
||||
namespace {
|
||||
|
||||
/* $OpenBSD: polevll.c,v 1.2 2013/11/12 20:35:09 martynas Exp $ */
|
||||
/*
|
||||
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
|
||||
*
|
||||
* Permission to use, copy, modify, and distribute this software for any
|
||||
* purpose with or without fee is hereby granted, provided that the above
|
||||
* copyright notice and this permission notice appear in all copies.
|
||||
*
|
||||
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
|
||||
* WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
|
||||
* MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
|
||||
* ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
|
||||
* WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
|
||||
* ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
|
||||
* OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
|
||||
*/
|
||||
|
||||
static long double P[8] = {
|
||||
4.212760487471622013093E-5L, 4.542931960608009155600E-4L,
|
||||
4.092666828394035500949E-3L, 2.385363243461108252554E-2L,
|
||||
1.113062816019361559013E-1L, 3.629515436640239168939E-1L,
|
||||
8.378004301573126728826E-1L, 1.000000000000000000009E0L,
|
||||
};
|
||||
static long double Q[9] = {
|
||||
-1.397148517476170440917E-5L, 2.346584059160635244282E-4L,
|
||||
-1.237799246653152231188E-3L, -7.955933682494738320586E-4L,
|
||||
2.773706565840072979165E-2L, -4.633887671244534213831E-2L,
|
||||
-2.243510905670329164562E-1L, 4.150160950588455434583E-1L,
|
||||
9.999999999999999999908E-1L,
|
||||
};
|
||||
|
||||
#define MAXGAML 1755.455L
|
||||
/*static const long double LOGPI = 1.14472988584940017414L;*/
|
||||
|
||||
/* Stirling's formula for the gamma function
|
||||
tgamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x))
|
||||
z(x) = x
|
||||
13 <= x <= 1024
|
||||
Relative error
|
||||
n=8, d=0
|
||||
Peak error = 9.44e-21
|
||||
Relative error spread = 8.8e-4
|
||||
*/
|
||||
|
||||
static long double STIR[9] = {
|
||||
7.147391378143610789273E-4L, -2.363848809501759061727E-5L,
|
||||
-5.950237554056330156018E-4L, 6.989332260623193171870E-5L,
|
||||
7.840334842744753003862E-4L, -2.294719747873185405699E-4L,
|
||||
-2.681327161876304418288E-3L, 3.472222222230075327854E-3L,
|
||||
8.333333333333331800504E-2L,
|
||||
};
|
||||
|
||||
#define MAXSTIR 1024.0L
|
||||
static const long double SQTPI = 2.50662827463100050242E0L;
|
||||
|
||||
/* 1/tgamma(x) = z P(z)
|
||||
* z(x) = 1/x
|
||||
* 0 < x < 0.03125
|
||||
* Peak relative error 4.2e-23
|
||||
*/
|
||||
|
||||
static long double S[9] = {
|
||||
-1.193945051381510095614E-3L, 7.220599478036909672331E-3L,
|
||||
-9.622023360406271645744E-3L, -4.219773360705915470089E-2L,
|
||||
1.665386113720805206758E-1L, -4.200263503403344054473E-2L,
|
||||
-6.558780715202540684668E-1L, 5.772156649015328608253E-1L,
|
||||
1.000000000000000000000E0L,
|
||||
};
|
||||
|
||||
/* 1/tgamma(-x) = z P(z)
|
||||
* z(x) = 1/x
|
||||
* 0 < x < 0.03125
|
||||
* Peak relative error 5.16e-23
|
||||
* Relative error spread = 2.5e-24
|
||||
*/
|
||||
|
||||
static long double SN[9] = {
|
||||
1.133374167243894382010E-3L, 7.220837261893170325704E-3L,
|
||||
9.621911155035976733706E-3L, -4.219773343731191721664E-2L,
|
||||
-1.665386113944413519335E-1L, -4.200263503402112910504E-2L,
|
||||
6.558780715202536547116E-1L, 5.772156649015328608727E-1L,
|
||||
-1.000000000000000000000E0L,
|
||||
};
|
||||
|
||||
static const long double PIL = 3.1415926535897932384626L;
|
||||
|
||||
static long double stirf(long double);
|
||||
|
||||
long double __polevll(long double x, void *PP, int n)
|
||||
{
|
||||
long double y;
|
||||
long double *P;
|
||||
|
||||
P = (long double *)PP;
|
||||
y = *P++;
|
||||
do
|
||||
{
|
||||
y = y * x + *P++;
|
||||
} while (--n);
|
||||
|
||||
return (y);
|
||||
}
|
||||
|
||||
/* Gamma function computed by Stirling's formula.
|
||||
*/
|
||||
static long double stirf(long double x)
|
||||
{
|
||||
long double y, w, v;
|
||||
|
||||
w = 1.0L / x;
|
||||
/* For large x, use rational coefficients from the analytical expansion. */
|
||||
if (x > 1024.0L)
|
||||
w = (((((6.97281375836585777429E-5L * w + 7.84039221720066627474E-4L)
|
||||
* w
|
||||
- 2.29472093621399176955E-4L)
|
||||
* w
|
||||
- 2.68132716049382716049E-3L)
|
||||
* w
|
||||
+ 3.47222222222222222222E-3L)
|
||||
* w
|
||||
+ 8.33333333333333333333E-2L)
|
||||
* w
|
||||
+ 1.0L;
|
||||
else
|
||||
w = 1.0L + w * __polevll(w, STIR, 8);
|
||||
y = expl(x);
|
||||
if (x > MAXSTIR)
|
||||
{ /* Avoid overflow in pow() */
|
||||
v = powl(x, 0.5L * x - 0.25L);
|
||||
y = v * (v / y);
|
||||
}
|
||||
else
|
||||
{
|
||||
y = powl(x, x - 0.5L) / y;
|
||||
}
|
||||
y = SQTPI * y * w;
|
||||
return (y);
|
||||
}
|
||||
|
||||
long double olm_tgammal(long double x)
|
||||
{
|
||||
long double p, q, z;
|
||||
int i;
|
||||
|
||||
if (isnan(x)) return (NAN);
|
||||
if (x == INFINITY) return (INFINITY);
|
||||
if (x == -INFINITY) return (x - x);
|
||||
if (x == 0.0L) return (1.0L / x);
|
||||
q = fabsl(x);
|
||||
|
||||
if (q > 13.0L)
|
||||
{
|
||||
int sign = 1;
|
||||
if (q > MAXGAML) goto goverf;
|
||||
if (x < 0.0L)
|
||||
{
|
||||
p = floorl(q);
|
||||
if (p == q) return (x - x) / (x - x);
|
||||
i = p;
|
||||
if ((i & 1) == 0) sign = -1;
|
||||
z = q - p;
|
||||
if (z > 0.5L)
|
||||
{
|
||||
p += 1.0L;
|
||||
z = q - p;
|
||||
}
|
||||
z = q * sinl(PIL * z);
|
||||
z = fabsl(z) * stirf(q);
|
||||
if (z <= PIL / LDBL_MAX)
|
||||
{
|
||||
goverf:
|
||||
return (sign * INFINITY);
|
||||
}
|
||||
z = PIL / z;
|
||||
}
|
||||
else
|
||||
{
|
||||
z = stirf(x);
|
||||
}
|
||||
return (sign * z);
|
||||
}
|
||||
|
||||
z = 1.0L;
|
||||
while (x >= 3.0L)
|
||||
{
|
||||
x -= 1.0L;
|
||||
z *= x;
|
||||
}
|
||||
|
||||
while (x < -0.03125L)
|
||||
{
|
||||
z /= x;
|
||||
x += 1.0L;
|
||||
}
|
||||
|
||||
if (x <= 0.03125L) goto smallarg;
|
||||
|
||||
while (x < 2.0L)
|
||||
{
|
||||
z /= x;
|
||||
x += 1.0L;
|
||||
}
|
||||
|
||||
if (x == 2.0L) return (z);
|
||||
|
||||
x -= 2.0L;
|
||||
p = __polevll(x, P, 7);
|
||||
q = __polevll(x, Q, 8);
|
||||
z = z * p / q;
|
||||
return z;
|
||||
|
||||
smallarg:
|
||||
if (x == 0.0L)
|
||||
return (x - x) / (x - x);
|
||||
else
|
||||
{
|
||||
if (x < 0.0L)
|
||||
{
|
||||
x = -x;
|
||||
q = z / (x * __polevll(x, SN, 8));
|
||||
}
|
||||
else
|
||||
q = z / (x * __polevll(x, S, 8));
|
||||
}
|
||||
return q;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
double reference_tgamma(double x)
|
||||
{
|
||||
if (x < 0.0f)
|
||||
{
|
||||
float fx = x;
|
||||
if (fx == std::floor(fx))
|
||||
{
|
||||
return std::numeric_limits<double>::signaling_NaN();
|
||||
}
|
||||
|
||||
// The gamma function satisfies the reflection formula:
|
||||
// tgamma(x) = pi / (sin(pi*x) * tgamma(1-x))
|
||||
// For large negative non-integer arguments, tgamma(1-x) should overflow
|
||||
// to +inf, causing the reflection formula pi/(sin(pi*x)*tgamma(1-x)) to
|
||||
// collapse to +/-0. However, reference olm_tgamma does not handle this
|
||||
// case correctly, so it must be handled externally. The sign is
|
||||
// determined analytically from the alternating sign rule: for x in
|
||||
// (-(n+1), -n), sign(tgamma(x)) = (-1)^(n+1). The threshold at which
|
||||
// tgamma(1-x) overflows float128 (~10^4932 max) is around argument
|
||||
// 1756, so x < -1755 is sufficient to trigger this path.
|
||||
|
||||
if (x < -1755.0f)
|
||||
{
|
||||
long n = (long)std::floor(-fx);
|
||||
bool negative = (n % 2 == 0);
|
||||
return negative ? -0.0 : 0.0;
|
||||
}
|
||||
}
|
||||
return olm_tgammal(x);
|
||||
}
|
||||
|
||||
#pragma mark -
|
||||
#pragma mark Double testing
|
||||
|
||||
|
|
@ -5815,3 +6077,32 @@ long double reference_erfl(long double x) { return erf(x); }
|
|||
|
||||
double reference_erfc(double x) { return erfc(x); }
|
||||
double reference_erf(double x) { return erf(x); }
|
||||
long double reference_tgammal(long double x)
|
||||
{
|
||||
if (x < 0.0)
|
||||
{
|
||||
if (x == std::floor(x))
|
||||
{
|
||||
return std::numeric_limits<long double>::signaling_NaN();
|
||||
}
|
||||
|
||||
// The gamma function satisfies the reflection formula:
|
||||
// tgamma(x) = pi / (sin(pi*x) * tgamma(1-x))
|
||||
// For large negative non-integer arguments, tgamma(1-x) should overflow
|
||||
// to +inf, causing the reflection formula pi/(sin(pi*x)*tgamma(1-x)) to
|
||||
// collapse to +/-0. However, reference olm_tgamma does not handle this
|
||||
// case correctly, so it must be handled externally. The sign is
|
||||
// determined analytically from the alternating sign rule: for x in
|
||||
// (-(n+1), -n), sign(tgamma(x)) = (-1)^(n+1). The threshold at which
|
||||
// tgamma(1-x) overflows float128 (~10^4932 max) is around argument
|
||||
// 1756, so x < -1755 is sufficient to trigger this path.
|
||||
|
||||
if (x < -1755.0)
|
||||
{
|
||||
long n = (long)std::floor(-x);
|
||||
bool negative = (n % 2 == 0);
|
||||
return negative ? -0.0L : 0.0L;
|
||||
}
|
||||
}
|
||||
return olm_tgammal(x);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -116,6 +116,11 @@ double reference_ldexp(double x, int n);
|
|||
|
||||
double reference_assignment(double x);
|
||||
int reference_not(double x);
|
||||
|
||||
double reference_erfc(double x);
|
||||
double reference_erf(double x);
|
||||
double reference_tgamma(double x);
|
||||
|
||||
// -- for testing fast-relaxed
|
||||
|
||||
double reference_relaxed_acos(double);
|
||||
|
|
@ -238,6 +243,6 @@ int reference_notl(long double x);
|
|||
|
||||
long double reference_erfcl(long double x);
|
||||
long double reference_erfl(long double x);
|
||||
double reference_erfc(double x);
|
||||
double reference_erf(double x);
|
||||
long double reference_tgammal(long double x);
|
||||
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -398,6 +398,8 @@ int TestFunc_Double_Double(const Func *f, MTdata d, bool relaxedMode)
|
|||
}
|
||||
}
|
||||
|
||||
test_info.isRangeLimited = 0;
|
||||
|
||||
// Init the kernels
|
||||
BuildKernelInfo build_info{ test_info.threadCount, test_info.k,
|
||||
test_info.programs, f->nameInCode,
|
||||
|
|
|
|||
|
|
@ -78,6 +78,7 @@ struct TestInfo
|
|||
int isRangeLimited; // 1 if the function is only to be evaluated over a
|
||||
// range
|
||||
float half_sin_cos_tan_limit;
|
||||
|
||||
bool relaxedMode; // True if test is running in relaxed mode, false
|
||||
// otherwise.
|
||||
};
|
||||
|
|
@ -352,14 +353,17 @@ cl_int Test(cl_uint job_id, cl_uint thread_id, void *data)
|
|||
}
|
||||
|
||||
// half_sin/cos/tan are only valid between +-2**16, Inf, NaN
|
||||
if (isRangeLimited
|
||||
&& fabsf(s[j]) > MAKE_HEX_FLOAT(0x1.0p16f, 0x1L, 16)
|
||||
&& fabsf(s[j]) < INFINITY)
|
||||
if (isRangeLimited)
|
||||
{
|
||||
if (fabsf(test) <= half_sin_cos_tan_limit)
|
||||
if (half_sin_cos_tan_limit > 0
|
||||
&& fabsf(s[j]) > MAKE_HEX_FLOAT(0x1.0p16f, 0x1L, 16)
|
||||
&& fabsf(s[j]) < INFINITY)
|
||||
{
|
||||
err = 0;
|
||||
fail = 0;
|
||||
if (fabsf(test) <= half_sin_cos_tan_limit)
|
||||
{
|
||||
err = 0;
|
||||
fail = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -78,6 +78,7 @@ cl_int TestHalf(cl_uint job_id, cl_uint thread_id, void *data)
|
|||
|
||||
int isRangeLimited = job->isRangeLimited;
|
||||
float half_sin_cos_tan_limit = job->half_sin_cos_tan_limit;
|
||||
|
||||
int ftz = job->ftz;
|
||||
|
||||
std::vector<float> s(0);
|
||||
|
|
@ -229,14 +230,18 @@ cl_int TestHalf(cl_uint job_id, cl_uint thread_id, void *data)
|
|||
int fail = !(fabsf(err) <= ulps);
|
||||
|
||||
// half_sin/cos/tan are only valid between +-2**16, Inf, NaN
|
||||
if (isRangeLimited
|
||||
&& fabsf(s[j]) > MAKE_HEX_FLOAT(0x1.0p16f, 0x1L, 16)
|
||||
&& fabsf(s[j]) < INFINITY)
|
||||
if (isRangeLimited)
|
||||
{
|
||||
if (fabsf(test) <= half_sin_cos_tan_limit)
|
||||
|
||||
if (half_sin_cos_tan_limit > 0
|
||||
&& fabsf(s[j]) > MAKE_HEX_FLOAT(0x1.0p16f, 0x1L, 16)
|
||||
&& fabsf(s[j]) < INFINITY)
|
||||
{
|
||||
err = 0;
|
||||
fail = 0;
|
||||
if (fabsf(test) <= half_sin_cos_tan_limit)
|
||||
{
|
||||
err = 0;
|
||||
fail = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
Loading…
Reference in a new issue