[Libclc-dev] [PATCH 2/3] remquo: Port from amd builtins
Jan Vesely via Libclc-dev
libclc-dev at lists.llvm.org
Tue May 1 15:01:26 PDT 2018
double version passes on carrizo. float version fails on denormals.
Signed-off-by: Jan Vesely <jan.vesely at rutgers.edu>
---
denormal failure is fixed in the next patch
generic/include/clc/clc.h | 1 +
generic/include/clc/math/remquo.h | 18 +++
generic/include/clc/math/remquo.inc | 1 +
generic/include/math/clc_remquo.h | 8 ++
generic/lib/SOURCES | 2 +
generic/lib/math/clc_remquo.cl | 254 ++++++++++++++++++++++++++++++++++++
generic/lib/math/remquo.cl | 17 +++
generic/lib/math/remquo.inc | 6 +
8 files changed, 307 insertions(+)
create mode 100644 generic/include/clc/math/remquo.h
create mode 100644 generic/include/clc/math/remquo.inc
create mode 100644 generic/include/math/clc_remquo.h
create mode 100644 generic/lib/math/clc_remquo.cl
create mode 100644 generic/lib/math/remquo.cl
create mode 100644 generic/lib/math/remquo.inc
diff --git a/generic/include/clc/clc.h b/generic/include/clc/clc.h
index 75ca7d0..171b06a 100644
--- a/generic/include/clc/clc.h
+++ b/generic/include/clc/clc.h
@@ -104,6 +104,7 @@
#include <clc/math/pown.h>
#include <clc/math/powr.h>
#include <clc/math/remainder.h>
+#include <clc/math/remquo.h>
#include <clc/math/rint.h>
#include <clc/math/rootn.h>
#include <clc/math/round.h>
diff --git a/generic/include/clc/math/remquo.h b/generic/include/clc/math/remquo.h
new file mode 100644
index 0000000..7daf82f
--- /dev/null
+++ b/generic/include/clc/math/remquo.h
@@ -0,0 +1,18 @@
+#define __CLC_FUNCTION remquo
+
+#define __CLC_BODY <clc/math/remquo.inc>
+#define __CLC_ADDRESS_SPACE global
+#include <clc/math/gentype.inc>
+#undef __CLC_ADDRESS_SPACE
+
+#define __CLC_BODY <clc/math/remquo.inc>
+#define __CLC_ADDRESS_SPACE local
+#include <clc/math/gentype.inc>
+#undef __CLC_ADDRESS_SPACE
+
+#define __CLC_BODY <clc/math/remquo.inc>
+#define __CLC_ADDRESS_SPACE private
+#include <clc/math/gentype.inc>
+#undef __CLC_ADDRESS_SPACE
+
+#undef __CLC_FUNCTION
diff --git a/generic/include/clc/math/remquo.inc b/generic/include/clc/math/remquo.inc
new file mode 100644
index 0000000..42c7b67
--- /dev/null
+++ b/generic/include/clc/math/remquo.inc
@@ -0,0 +1 @@
+_CLC_OVERLOAD _CLC_DECL __CLC_GENTYPE __CLC_FUNCTION(__CLC_GENTYPE x, __CLC_GENTYPE y, __CLC_ADDRESS_SPACE __CLC_INTN *q);
diff --git a/generic/include/math/clc_remquo.h b/generic/include/math/clc_remquo.h
new file mode 100644
index 0000000..ed57ec9
--- /dev/null
+++ b/generic/include/math/clc_remquo.h
@@ -0,0 +1,8 @@
+#define __CLC_FUNCTION __clc_remquo
+
+#define __CLC_BODY <clc/math/remquo.inc>
+#define __CLC_ADDRESS_SPACE private
+#include <clc/math/gentype.inc>
+#undef __CLC_ADDRESS_SPACE
+
+#undef __CLC_FUNCTION
diff --git a/generic/lib/SOURCES b/generic/lib/SOURCES
index 8598fa5..1f1a427 100644
--- a/generic/lib/SOURCES
+++ b/generic/lib/SOURCES
@@ -163,6 +163,8 @@ math/clc_powr.cl
math/powr.cl
math/clc_remainder.cl
math/remainder.cl
+math/clc_remquo.cl
+math/remquo.cl
math/clc_rootn.cl
math/rootn.cl
math/sin.cl
diff --git a/generic/lib/math/clc_remquo.cl b/generic/lib/math/clc_remquo.cl
new file mode 100644
index 0000000..ff741f6
--- /dev/null
+++ b/generic/lib/math/clc_remquo.cl
@@ -0,0 +1,254 @@
+/*
+ * Copyright (c) 2014 Advanced Micro Devices, Inc.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ */
+
+#include <clc/clc.h>
+
+#include <math/clc_remainder.h>
+#include "../clcmacro.h"
+#include "config.h"
+#include "math.h"
+
+_CLC_DEF _CLC_OVERLOAD float __clc_remquo(float x, float y, __private int *quo)
+{
+ int ux = as_int(x);
+ int ax = ux & EXSIGNBIT_SP32;
+ float xa = as_float(ax);
+ int sx = ux ^ ax;
+ int ex = ax >> EXPSHIFTBITS_SP32;
+
+ int uy = as_int(y);
+ int ay = uy & EXSIGNBIT_SP32;
+ float ya = as_float(ay);
+ int sy = uy ^ ay;
+ int ey = ay >> EXPSHIFTBITS_SP32;
+
+ float xr = as_float(0x3f800000 | (ax & 0x007fffff));
+ float yr = as_float(0x3f800000 | (ay & 0x007fffff));
+ int c;
+ int k = ex - ey;
+
+ uint q = 0;
+
+ while (k > 0) {
+ c = xr >= yr;
+ q = (q << 1) | c;
+ xr -= c ? yr : 0.0f;
+ xr += xr;
+ --k;
+ }
+
+ c = xr > yr;
+ q = (q << 1) | c;
+ xr -= c ? yr : 0.0f;
+
+ int lt = ex < ey;
+
+ q = lt ? 0 : q;
+ xr = lt ? xa : xr;
+ yr = lt ? ya : yr;
+
+ c = (yr < 2.0f * xr) | ((yr == 2.0f * xr) & ((q & 0x1) == 0x1));
+ xr -= c ? yr : 0.0f;
+ q += c;
+
+ float s = as_float(ey << EXPSHIFTBITS_SP32);
+ xr *= lt ? 1.0f : s;
+
+ int qsgn = sx == sy ? 1 : -1;
+ int quot = (q & 0x7f) * qsgn;
+
+ c = ax == ay;
+ quot = c ? qsgn : quot;
+ xr = c ? 0.0f : xr;
+
+ xr = as_float(sx ^ as_int(xr));
+
+ c = ax > PINFBITPATT_SP32 | ay > PINFBITPATT_SP32 | ax == PINFBITPATT_SP32 | ay == 0;
+ quot = c ? 0 : quot;
+ xr = c ? as_float(QNANBITPATT_SP32) : xr;
+
+ *quo = quot;
+
+ return xr;
+}
+// remquo singature is special, we don't have macro for this
+#define __VEC_REMQUO(TYPE, VEC_SIZE, HALF_VEC_SIZE) \
+_CLC_DEF _CLC_OVERLOAD TYPE##VEC_SIZE __clc_remquo(TYPE##VEC_SIZE x, TYPE##VEC_SIZE y, __private int##VEC_SIZE *quo) \
+{ \
+ int##HALF_VEC_SIZE lo, hi; \
+ TYPE##VEC_SIZE ret; \
+ ret.lo = __clc_remquo(x.lo, y.lo, &lo); \
+ ret.hi = __clc_remquo(x.hi, y.hi, &hi); \
+ (*quo).lo = lo; \
+ (*quo).hi = hi; \
+ return ret; \
+}
+__VEC_REMQUO(float, 2,)
+__VEC_REMQUO(float, 3, 2)
+__VEC_REMQUO(float, 4, 2)
+__VEC_REMQUO(float, 8, 4)
+__VEC_REMQUO(float, 16, 8)
+
+#ifdef cl_khr_fp64
+_CLC_DEF _CLC_OVERLOAD double __clc_remquo(double x, double y, __private int *pquo)
+{
+ ulong ux = as_ulong(x);
+ ulong ax = ux & ~SIGNBIT_DP64;
+ ulong xsgn = ux ^ ax;
+ double dx = as_double(ax);
+ int xexp = convert_int(ax >> EXPSHIFTBITS_DP64);
+ int xexp1 = 11 - (int) clz(ax & MANTBITS_DP64);
+ xexp1 = xexp < 1 ? xexp1 : xexp;
+
+ ulong uy = as_ulong(y);
+ ulong ay = uy & ~SIGNBIT_DP64;
+ double dy = as_double(ay);
+ int yexp = convert_int(ay >> EXPSHIFTBITS_DP64);
+ int yexp1 = 11 - (int) clz(ay & MANTBITS_DP64);
+ yexp1 = yexp < 1 ? yexp1 : yexp;
+
+ int qsgn = ((ux ^ uy) & SIGNBIT_DP64) == 0UL ? 1 : -1;
+
+ // First assume |x| > |y|
+
+ // Set ntimes to the number of times we need to do a
+ // partial remainder. If the exponent of x is an exact multiple
+ // of 53 larger than the exponent of y, and the mantissa of x is
+ // less than the mantissa of y, ntimes will be one too large
+ // but it doesn't matter - it just means that we'll go round
+ // the loop below one extra time.
+ int ntimes = max(0, (xexp1 - yexp1) / 53);
+ double w = ldexp(dy, ntimes * 53);
+ w = ntimes == 0 ? dy : w;
+ double scale = ntimes == 0 ? 1.0 : 0x1.0p-53;
+
+ // Each time round the loop we compute a partial remainder.
+ // This is done by subtracting a large multiple of w
+ // from x each time, where w is a scaled up version of y.
+ // The subtraction must be performed exactly in quad
+ // precision, though the result at each stage can
+ // fit exactly in a double precision number.
+ int i;
+ double t, v, p, pp;
+
+ for (i = 0; i < ntimes; i++) {
+ // Compute integral multiplier
+ t = trunc(dx / w);
+
+ // Compute w * t in quad precision
+ p = w * t;
+ pp = fma(w, t, -p);
+
+ // Subtract w * t from dx
+ v = dx - p;
+ dx = v + (((dx - v) - p) - pp);
+
+ // If t was one too large, dx will be negative. Add back one w.
+ dx += dx < 0.0 ? w : 0.0;
+
+ // Scale w down by 2^(-53) for the next iteration
+ w *= scale;
+ }
+
+ // One more time
+ // Variable todd says whether the integer t is odd or not
+ t = floor(dx / w);
+ long lt = (long)t;
+ int todd = lt & 1;
+
+ p = w * t;
+ pp = fma(w, t, -p);
+ v = dx - p;
+ dx = v + (((dx - v) - p) - pp);
+ i = dx < 0.0;
+ todd ^= i;
+ dx += i ? w : 0.0;
+
+ lt -= i;
+
+ // At this point, dx lies in the range [0,dy)
+
+ // For the remainder function, we need to adjust dx
+ // so that it lies in the range (-y/2, y/2] by carefully
+ // subtracting w (== dy == y) if necessary. The rigmarole
+ // with todd is to get the correct sign of the result
+ // when x/y lies exactly half way between two integers,
+ // when we need to choose the even integer.
+
+ int al = (2.0*dx > w) | (todd & (2.0*dx == w));
+ double dxl = dx - (al ? w : 0.0);
+
+ int ag = (dx > 0.5*w) | (todd & (dx == 0.5*w));
+ double dxg = dx - (ag ? w : 0.0);
+
+ dx = dy < 0x1.0p+1022 ? dxl : dxg;
+ lt += dy < 0x1.0p+1022 ? al : ag;
+ int quo = ((int)lt & 0x7f) * qsgn;
+
+ double ret = as_double(xsgn ^ as_ulong(dx));
+ dx = as_double(ax);
+
+ // Now handle |x| == |y|
+ int c = dx == dy;
+ t = as_double(xsgn);
+ quo = c ? qsgn : quo;
+ ret = c ? t : ret;
+
+ // Next, handle |x| < |y|
+ c = dx < dy;
+ quo = c ? 0 : quo;
+ ret = c ? x : ret;
+
+ c &= (yexp < 1023 & 2.0*dx > dy) | (dx > 0.5*dy);
+ quo = c ? qsgn : quo;
+ // we could use a conversion here instead since qsgn = +-1
+ p = qsgn == 1 ? -1.0 : 1.0;
+ t = fma(y, p, x);
+ ret = c ? t : ret;
+
+ // We don't need anything special for |x| == 0
+
+ // |y| is 0
+ c = dy == 0.0;
+ quo = c ? 0 : quo;
+ ret = c ? as_double(QNANBITPATT_DP64) : ret;
+
+ // y is +-Inf, NaN
+ c = yexp > BIASEDEMAX_DP64;
+ quo = c ? 0 : quo;
+ t = y == y ? x : y;
+ ret = c ? t : ret;
+
+ // x is +=Inf, NaN
+ c = xexp > BIASEDEMAX_DP64;
+ quo = c ? 0 : quo;
+ ret = c ? as_double(QNANBITPATT_DP64) : ret;
+
+ *pquo = quo;
+ return ret;
+}
+__VEC_REMQUO(double, 2,)
+__VEC_REMQUO(double, 3, 2)
+__VEC_REMQUO(double, 4, 2)
+__VEC_REMQUO(double, 8, 4)
+__VEC_REMQUO(double, 16, 8)
+#endif
diff --git a/generic/lib/math/remquo.cl b/generic/lib/math/remquo.cl
new file mode 100644
index 0000000..fc29b36
--- /dev/null
+++ b/generic/lib/math/remquo.cl
@@ -0,0 +1,17 @@
+#include <clc/clc.h>
+#include <math/clc_remquo.h>
+
+#define __CLC_BODY <remquo.inc>
+#define __CLC_ADDRESS_SPACE global
+#include <clc/math/gentype.inc>
+#undef __CLC_ADDRESS_SPACE
+
+#define __CLC_BODY <remquo.inc>
+#define __CLC_ADDRESS_SPACE local
+#include <clc/math/gentype.inc>
+#undef __CLC_ADDRESS_SPACE
+
+#define __CLC_BODY <remquo.inc>
+#define __CLC_ADDRESS_SPACE private
+#include <clc/math/gentype.inc>
+#undef __CLC_ADDRESS_SPACE
diff --git a/generic/lib/math/remquo.inc b/generic/lib/math/remquo.inc
new file mode 100644
index 0000000..c1de78a
--- /dev/null
+++ b/generic/lib/math/remquo.inc
@@ -0,0 +1,6 @@
+_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE remquo(__CLC_GENTYPE x, __CLC_GENTYPE y, __CLC_ADDRESS_SPACE __CLC_INTN *q) {
+ __CLC_INTN local_q;
+ __CLC_GENTYPE ret = __clc_remquo(x, y, &local_q);
+ *q = local_q;
+ return ret;
+}
--
2.14.3
More information about the Libclc-dev
mailing list