add built-in function "lgamma", "lgamma_r"

also include test cases

Signed-off-by: Homer Hsing <homer.xing@intel.com>
Reviewed-by: Zhigang Gong <zhigang.gong@linux.intel.com>

diff --git a/backend/src/builtin_vector_proto.def b/backend/src/builtin_vector_proto.def
index 2a057bb..2a3daf2 100644 (file)
--- 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 (file)
--- 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 (file)
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 (file)
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 6016938..69c7976 100644 (file)
--- 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 (file)
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 (file)
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);