[Beignet] [PATCH] add built-in function "lgamma", "lgamma_r"
Homer Hsing
homer.xing at intel.com
Sun Aug 25 21:51:53 PDT 2013
also include test cases
Signed-off-by: Homer Hsing <homer.xing at intel.com>
---
backend/src/builtin_vector_proto.def | 14 +-
backend/src/ocl_stdlib.tmpl.h | 384 +++++++++++++++++++++++++++++++++++
kernels/builtin_lgamma.cl | 4 +
kernels/builtin_lgamma_r.cl | 4 +
utests/CMakeLists.txt | 2 +
utests/builtin_lgamma.cpp | 40 ++++
utests/builtin_lgamma_r.cpp | 46 +++++
7 files changed, 487 insertions(+), 7 deletions(-)
create mode 100644 kernels/builtin_lgamma.cl
create mode 100644 kernels/builtin_lgamma_r.cl
create mode 100644 utests/builtin_lgamma.cpp
create mode 100644 utests/builtin_lgamma_r.cpp
diff --git a/backend/src/builtin_vector_proto.def b/backend/src/builtin_vector_proto.def
index 2a057bb..2a3daf2 100644
--- a/backend/src/builtin_vector_proto.def
+++ b/backend/src/builtin_vector_proto.def
@@ -61,13 +61,13 @@ float ldexp (float x, int k)
doublen ldexp (doublen x, intn k)
doublen ldexp (doublen x, int k)
double ldexp (double x, int k)
-#gentype lgamma (gentype x)
-#floatn lgamma_r (floatn x, __global intn *signp)
-#floatn lgamma_r (floatn x, __local intn *signp)
-#floatn lgamma_r (floatn x, __private intn *signp)
-#float lgamma_r (float x, __global int *signp)
-#float lgamma_r (float x, __local int *signp)
-#float lgamma_r (float x, __private int *signp)
+gentype lgamma (gentype x)
+floatn lgamma_r (floatn x, __global intn *signp)
+floatn lgamma_r (floatn x, __local intn *signp)
+floatn lgamma_r (floatn x, __private intn *signp)
+float lgamma_r (float x, __global int *signp)
+float lgamma_r (float x, __local int *signp)
+float lgamma_r (float x, __private int *signp)
#doublen lgamma_r (doublen x, __global intn *signp)
#doublen lgamma_r (doublen x, __local intn *signp)
#doublen lgamma_r (doublen x, __private intn *signp)
diff --git a/backend/src/ocl_stdlib.tmpl.h b/backend/src/ocl_stdlib.tmpl.h
index c8d20b6..ac1999d 100644
--- a/backend/src/ocl_stdlib.tmpl.h
+++ b/backend/src/ocl_stdlib.tmpl.h
@@ -823,6 +823,390 @@ INLINE_OVERLOADABLE float tgamma(float x) {
r = nadj - r;
return r;
}
+
+INLINE_OVERLOADABLE float lgamma(float x) {
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+ const float
+ zero= 0.,
+ one = 1.0000000000e+00,
+ pi = 3.1415927410e+00,
+ a0 = 7.7215664089e-02,
+ a1 = 3.2246702909e-01,
+ a2 = 6.7352302372e-02,
+ a3 = 2.0580807701e-02,
+ a4 = 7.3855509982e-03,
+ a5 = 2.8905137442e-03,
+ a6 = 1.1927076848e-03,
+ a7 = 5.1006977446e-04,
+ a8 = 2.2086278477e-04,
+ a9 = 1.0801156895e-04,
+ a10 = 2.5214456400e-05,
+ a11 = 4.4864096708e-05,
+ tc = 1.4616321325e+00,
+ tf = -1.2148628384e-01,
+ tt = 6.6971006518e-09,
+ t0 = 4.8383611441e-01,
+ t1 = -1.4758771658e-01,
+ t2 = 6.4624942839e-02,
+ t3 = -3.2788541168e-02,
+ t4 = 1.7970675603e-02,
+ t5 = -1.0314224288e-02,
+ t6 = 6.1005386524e-03,
+ t7 = -3.6845202558e-03,
+ t8 = 2.2596477065e-03,
+ t9 = -1.4034647029e-03,
+ t10 = 8.8108185446e-04,
+ t11 = -5.3859531181e-04,
+ t12 = 3.1563205994e-04,
+ t13 = -3.1275415677e-04,
+ t14 = 3.3552918467e-04,
+ u0 = -7.7215664089e-02,
+ u1 = 6.3282704353e-01,
+ u2 = 1.4549225569e+00,
+ u3 = 9.7771751881e-01,
+ u4 = 2.2896373272e-01,
+ u5 = 1.3381091878e-02,
+ v1 = 2.4559779167e+00,
+ v2 = 2.1284897327e+00,
+ v3 = 7.6928514242e-01,
+ v4 = 1.0422264785e-01,
+ v5 = 3.2170924824e-03,
+ s0 = -7.7215664089e-02,
+ s1 = 2.1498242021e-01,
+ s2 = 3.2577878237e-01,
+ s3 = 1.4635047317e-01,
+ s4 = 2.6642270386e-02,
+ s5 = 1.8402845599e-03,
+ s6 = 3.1947532989e-05,
+ r1 = 1.3920053244e+00,
+ r2 = 7.2193557024e-01,
+ r3 = 1.7193385959e-01,
+ r4 = 1.8645919859e-02,
+ r5 = 7.7794247773e-04,
+ r6 = 7.3266842264e-06,
+ w0 = 4.1893854737e-01,
+ w1 = 8.3333335817e-02,
+ w2 = -2.7777778450e-03,
+ w3 = 7.9365057172e-04,
+ w4 = -5.9518753551e-04,
+ w5 = 8.3633989561e-04,
+ w6 = -1.6309292987e-03;
+ float t, y, z, nadj, p, p1, p2, p3, q, r, w;
+ int i, hx, ix;
+ nadj = 0;
+ hx = *(int *)&x;
+ ix = hx & 0x7fffffff;
+ if (ix >= 0x7f800000)
+ return x * x;
+ if (ix == 0)
+ return ((x + one) / zero);
+ if (ix < 0x1c800000) {
+ if (hx < 0) {
+ return -native_log(-x);
+ } else
+ return -native_log(x);
+ }
+ if (hx < 0) {
+ if (ix >= 0x4b000000)
+ return ((-x) / zero);
+ t = __gen_ocl_internal_sinpi(x);
+ if (t == zero)
+ return ((-x) / zero);
+ nadj = native_log(pi / __gen_ocl_fabs(t * x));
+ x = -x;
+ }
+ if (ix == 0x3f800000 || ix == 0x40000000)
+ r = 0;
+ else if (ix < 0x40000000) {
+ if (ix <= 0x3f666666) {
+ r = -native_log(x);
+ if (ix >= 0x3f3b4a20) {
+ y = one - x;
+ i = 0;
+ } else if (ix >= 0x3e6d3308) {
+ y = x - (tc - one);
+ i = 1;
+ } else {
+ y = x;
+ i = 2;
+ }
+ } else {
+ r = zero;
+ if (ix >= 0x3fdda618) {
+ y = (float) 2.0 - x;
+ i = 0;
+ }
+ else if (ix >= 0x3F9da620) {
+ y = x - tc;
+ i = 1;
+ }
+ else {
+ y = x - one;
+ i = 2;
+ }
+ }
+ switch (i) {
+ case 0:
+ z = y * y;
+ p1 = a0 + z * (a2 + z * (a4 + z * (a6 + z * (a8 + z * a10))));
+ p2 = z * (a1 + z * (a3 + z * (a5 + z * (a7 + z * (a9 + z * a11)))));
+ p = y * p1 + p2;
+ r += (p - (float) 0.5 * y);
+ break;
+ case 1:
+ z = y * y;
+ w = z * y;
+ p1 = t0 + w * (t3 + w * (t6 + w * (t9 + w * t12)));
+ p2 = t1 + w * (t4 + w * (t7 + w * (t10 + w * t13)));
+ p3 = t2 + w * (t5 + w * (t8 + w * (t11 + w * t14)));
+ p = z * p1 - (tt - w * (p2 + y * p3));
+ r += (tf + p);
+ break;
+ case 2:
+ p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * u5)))));
+ p2 = one + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * v5))));
+ r += (-(float) 0.5 * y + p1 / p2);
+ }
+ } else if (ix < 0x41000000) {
+ i = (int) x;
+ t = zero;
+ y = x - (float) i;
+ p = y * (s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))))));
+ q = one + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6)))));
+ r = .5f * y + p / q;
+ z = one;
+ switch (i) {
+ case 7:
+ z *= (y + (float) 6.0);
+ case 6:
+ z *= (y + (float) 5.0);
+ case 5:
+ z *= (y + (float) 4.0);
+ case 4:
+ z *= (y + (float) 3.0);
+ case 3:
+ z *= (y + (float) 2.0);
+ r += native_log(z);
+ break;
+ }
+
+ } else if (ix < 0x5c800000) {
+ t = native_log(x);
+ z = one / x;
+ y = z * z;
+ w = w0 + z * (w1 + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * w6)))));
+ r = (x - .5f) * (t - one) + w;
+ } else
+ r = x * (native_log(x) - one);
+ if (hx < 0)
+ r = nadj - r;
+ return r;
+}
+
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+#define BODY \
+ const float \
+ zero= 0., \
+ one = 1.0000000000e+00, \
+ pi = 3.1415927410e+00, \
+ a0 = 7.7215664089e-02, \
+ a1 = 3.2246702909e-01, \
+ a2 = 6.7352302372e-02, \
+ a3 = 2.0580807701e-02, \
+ a4 = 7.3855509982e-03, \
+ a5 = 2.8905137442e-03, \
+ a6 = 1.1927076848e-03, \
+ a7 = 5.1006977446e-04, \
+ a8 = 2.2086278477e-04, \
+ a9 = 1.0801156895e-04, \
+ a10 = 2.5214456400e-05, \
+ a11 = 4.4864096708e-05, \
+ tc = 1.4616321325e+00, \
+ tf = -1.2148628384e-01, \
+ tt = 6.6971006518e-09, \
+ t0 = 4.8383611441e-01, \
+ t1 = -1.4758771658e-01, \
+ t2 = 6.4624942839e-02, \
+ t3 = -3.2788541168e-02, \
+ t4 = 1.7970675603e-02, \
+ t5 = -1.0314224288e-02, \
+ t6 = 6.1005386524e-03, \
+ t7 = -3.6845202558e-03, \
+ t8 = 2.2596477065e-03, \
+ t9 = -1.4034647029e-03, \
+ t10 = 8.8108185446e-04, \
+ t11 = -5.3859531181e-04, \
+ t12 = 3.1563205994e-04, \
+ t13 = -3.1275415677e-04, \
+ t14 = 3.3552918467e-04, \
+ u0 = -7.7215664089e-02, \
+ u1 = 6.3282704353e-01, \
+ u2 = 1.4549225569e+00, \
+ u3 = 9.7771751881e-01, \
+ u4 = 2.2896373272e-01, \
+ u5 = 1.3381091878e-02, \
+ v1 = 2.4559779167e+00, \
+ v2 = 2.1284897327e+00, \
+ v3 = 7.6928514242e-01, \
+ v4 = 1.0422264785e-01, \
+ v5 = 3.2170924824e-03, \
+ s0 = -7.7215664089e-02, \
+ s1 = 2.1498242021e-01, \
+ s2 = 3.2577878237e-01, \
+ s3 = 1.4635047317e-01, \
+ s4 = 2.6642270386e-02, \
+ s5 = 1.8402845599e-03, \
+ s6 = 3.1947532989e-05, \
+ r1 = 1.3920053244e+00, \
+ r2 = 7.2193557024e-01, \
+ r3 = 1.7193385959e-01, \
+ r4 = 1.8645919859e-02, \
+ r5 = 7.7794247773e-04, \
+ r6 = 7.3266842264e-06, \
+ w0 = 4.1893854737e-01, \
+ w1 = 8.3333335817e-02, \
+ w2 = -2.7777778450e-03, \
+ w3 = 7.9365057172e-04, \
+ w4 = -5.9518753551e-04, \
+ w5 = 8.3633989561e-04, \
+ w6 = -1.6309292987e-03; \
+ float t, y, z, nadj, p, p1, p2, p3, q, r, w; \
+ int i, hx, ix; \
+ nadj = 0; \
+ hx = *(int *)&x; \
+ *signgamp = 1; \
+ ix = hx & 0x7fffffff; \
+ if (ix >= 0x7f800000) \
+ return x * x; \
+ if (ix == 0) \
+ return ((x + one) / zero); \
+ if (ix < 0x1c800000) { \
+ if (hx < 0) { \
+ *signgamp = -1; \
+ return -native_log(-x); \
+ } else \
+ return -native_log(x); \
+ } \
+ if (hx < 0) { \
+ if (ix >= 0x4b000000) \
+ return ((-x) / zero); \
+ t = __gen_ocl_internal_sinpi(x); \
+ if (t == zero) \
+ return ((-x) / zero); \
+ nadj = native_log(pi / __gen_ocl_fabs(t * x)); \
+ if (t < zero) \
+ *signgamp = -1; \
+ x = -x; \
+ } \
+ if (ix == 0x3f800000 || ix == 0x40000000) \
+ r = 0; \
+ else if (ix < 0x40000000) { \
+ if (ix <= 0x3f666666) { \
+ r = -native_log(x); \
+ if (ix >= 0x3f3b4a20) { \
+ y = one - x; \
+ i = 0; \
+ } else if (ix >= 0x3e6d3308) { \
+ y = x - (tc - one); \
+ i = 1; \
+ } else { \
+ y = x; \
+ i = 2; \
+ } \
+ } else { \
+ r = zero; \
+ if (ix >= 0x3fdda618) { \
+ y = (float) 2.0 - x; \
+ i = 0; \
+ } \
+ else if (ix >= 0x3F9da620) { \
+ y = x - tc; \
+ i = 1; \
+ } \
+ else { \
+ y = x - one; \
+ i = 2; \
+ } \
+ } \
+ switch (i) { \
+ case 0: \
+ z = y * y; \
+ p1 = a0 + z * (a2 + z * (a4 + z * (a6 + z * (a8 + z * a10)))); \
+ p2 = z * (a1 + z * (a3 + z * (a5 + z * (a7 + z * (a9 + z * a11))))); \
+ p = y * p1 + p2; \
+ r += (p - (float) 0.5 * y); \
+ break; \
+ case 1: \
+ z = y * y; \
+ w = z * y; \
+ p1 = t0 + w * (t3 + w * (t6 + w * (t9 + w * t12))); \
+ p2 = t1 + w * (t4 + w * (t7 + w * (t10 + w * t13))); \
+ p3 = t2 + w * (t5 + w * (t8 + w * (t11 + w * t14))); \
+ p = z * p1 - (tt - w * (p2 + y * p3)); \
+ r += (tf + p); \
+ break; \
+ case 2: \
+ p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * u5))))); \
+ p2 = one + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * v5)))); \
+ r += (-(float) 0.5 * y + p1 / p2); \
+ } \
+ } else if (ix < 0x41000000) { \
+ i = (int) x; \
+ t = zero; \
+ y = x - (float) i; \
+ p = y * (s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6)))))); \
+ q = one + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6))))); \
+ r = .5f * y + p / q; \
+ z = one; \
+ switch (i) { \
+ case 7: \
+ z *= (y + (float) 6.0); \
+ case 6: \
+ z *= (y + (float) 5.0); \
+ case 5: \
+ z *= (y + (float) 4.0); \
+ case 4: \
+ z *= (y + (float) 3.0); \
+ case 3: \
+ z *= (y + (float) 2.0); \
+ r += native_log(z); \
+ break; \
+ } \
+ \
+ } else if (ix < 0x5c800000) { \
+ t = native_log(x); \
+ z = one / x; \
+ y = z * z; \
+ w = w0 + z * (w1 + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * w6))))); \
+ r = (x - .5f) * (t - one) + w; \
+ } else \
+ r = x * (native_log(x) - one); \
+ if (hx < 0) \
+ r = nadj - r; \
+ return r;
+INLINE_OVERLOADABLE float lgamma_r(float x, global int *signgamp) { BODY; }
+INLINE_OVERLOADABLE float lgamma_r(float x, local int *signgamp) { BODY; }
+INLINE_OVERLOADABLE float lgamma_r(float x, private int *signgamp) { BODY; }
+#undef BODY
+
INLINE_OVERLOADABLE float native_log10(float x) {
return native_log2(x) * 0.3010299956f;
}
diff --git a/kernels/builtin_lgamma.cl b/kernels/builtin_lgamma.cl
new file mode 100644
index 0000000..85bf859
--- /dev/null
+++ b/kernels/builtin_lgamma.cl
@@ -0,0 +1,4 @@
+kernel void builtin_lgamma(global float *src, global float *dst) {
+ int i = get_global_id(0);
+ dst[i] = lgamma(src[i]);
+};
diff --git a/kernels/builtin_lgamma_r.cl b/kernels/builtin_lgamma_r.cl
new file mode 100644
index 0000000..71fcc36
--- /dev/null
+++ b/kernels/builtin_lgamma_r.cl
@@ -0,0 +1,4 @@
+kernel void builtin_lgamma_r(global float *src, global float *dst, global int *signp) {
+ int i = get_global_id(0);
+ dst[i] = lgamma_r(src[i], signp+i);
+};
diff --git a/utests/CMakeLists.txt b/utests/CMakeLists.txt
index 949c6b7..08b4e32 100644
--- a/utests/CMakeLists.txt
+++ b/utests/CMakeLists.txt
@@ -114,6 +114,8 @@ set (utests_sources
builtin_shuffle2.cpp
builtin_sign.cpp
builtin_sinpi.cpp
+ builtin_lgamma.cpp
+ builtin_lgamma_r.cpp
builtin_tgamma.cpp
buildin_work_dim.cpp
builtin_global_size.cpp
diff --git a/utests/builtin_lgamma.cpp b/utests/builtin_lgamma.cpp
new file mode 100644
index 0000000..876699a
--- /dev/null
+++ b/utests/builtin_lgamma.cpp
@@ -0,0 +1,40 @@
+#include <cmath>
+#include "utest_helper.hpp"
+
+void builtin_lgamma(void) {
+ const int n = 1024;
+ float src[n];
+
+ // Setup kernel and buffers
+ OCL_CREATE_KERNEL("builtin_lgamma");
+ OCL_CREATE_BUFFER(buf[0], 0, n * sizeof(float), NULL);
+ OCL_CREATE_BUFFER(buf[1], 0, n * sizeof(float), NULL);
+ OCL_SET_ARG(0, sizeof(cl_mem), &buf[0]);
+ OCL_SET_ARG(1, sizeof(cl_mem), &buf[1]);
+ globals[0] = n;
+ locals[0] = 16;
+
+ for (int j = 0; j < 1024; j++) {
+ OCL_MAP_BUFFER(0);
+ for (int i = 0; i < n; ++i) {
+ src[i] = ((float*) buf_data[0])[i] = (j * n + i + 1) * 0.001f;
+ }
+ OCL_UNMAP_BUFFER(0);
+
+ OCL_NDRANGE(1);
+
+ OCL_MAP_BUFFER(1);
+ float *dst = (float*) buf_data[1];
+ for (int i = 0; i < n; ++i) {
+ float cpu = lgamma(src[i]);
+ float gpu = dst[i];
+ if (fabsf(cpu - gpu) >= 1e-3) {
+ printf("%f %f %f\n", src[i], cpu, gpu);
+ OCL_ASSERT(0);
+ }
+ }
+ OCL_UNMAP_BUFFER(1);
+ }
+}
+
+MAKE_UTEST_FROM_FUNCTION (builtin_lgamma);
diff --git a/utests/builtin_lgamma_r.cpp b/utests/builtin_lgamma_r.cpp
new file mode 100644
index 0000000..b6e5d0e
--- /dev/null
+++ b/utests/builtin_lgamma_r.cpp
@@ -0,0 +1,46 @@
+#include <cmath>
+#include "utest_helper.hpp"
+
+void builtin_lgamma_r(void) {
+ const int n = 1024;
+ float src[n];
+
+ // Setup kernel and buffers
+ OCL_CREATE_KERNEL("builtin_lgamma_r");
+ OCL_CREATE_BUFFER(buf[0], 0, n * sizeof(float), NULL);
+ OCL_CREATE_BUFFER(buf[1], 0, n * sizeof(float), NULL);
+ OCL_CREATE_BUFFER(buf[2], 0, n * sizeof(int), NULL);
+ OCL_SET_ARG(0, sizeof(cl_mem), &buf[0]);
+ OCL_SET_ARG(1, sizeof(cl_mem), &buf[1]);
+ OCL_SET_ARG(2, sizeof(cl_mem), &buf[2]);
+ globals[0] = n;
+ locals[0] = 16;
+
+ for (int j = 0; j < 1024; j++) {
+ OCL_MAP_BUFFER(0);
+ for (int i = 0; i < n; ++i) {
+ src[i] = ((float*) buf_data[0])[i] = (j * n + i + 1) * 0.001f;
+ }
+ OCL_UNMAP_BUFFER(0);
+
+ OCL_NDRANGE(1);
+
+ OCL_MAP_BUFFER(1);
+ OCL_MAP_BUFFER(2);
+ float *dst = (float*) buf_data[1];
+ for (int i = 0; i < n; ++i) {
+ int cpu_signp;
+ float cpu = lgamma_r(src[i], &cpu_signp);
+ int gpu_signp = ((int*)buf_data[2])[i];
+ float gpu = dst[i];
+ if (cpu_signp != gpu_signp || fabsf(cpu - gpu) >= 1e-3) {
+ printf("%f %f %f\n", src[i], cpu, gpu);
+ OCL_ASSERT(0);
+ }
+ }
+ OCL_UNMAP_BUFFER(1);
+ OCL_UNMAP_BUFFER(2);
+ }
+}
+
+MAKE_UTEST_FROM_FUNCTION (builtin_lgamma_r);
--
1.8.1.2
More information about the Beignet
mailing list