[Libclc-dev] [PATCH] math: Implement remainder(x, y)
Aaron Watry via Libclc-dev
libclc-dev at lists.llvm.org
Wed Jan 18 20:04:29 PST 2017
Mostly ported from the amd-builtins branch.
The amd-builtins branch uses __amdil_improved_fdiv_f32 and FTZ which aren't available in generic CLC.
__amdil_improved_fdiv_f32 points to native_divide which does native_recip(y)*x.
Since we don't have native_divide or native_recip yet, I've just stuck an actual division here.
I've taken a shot at a replacement for FTZ(x), but feel free to suggest alternatives.
Tested via piglit on a Radeon HD 7850 using the tests just sent to that list.
v2: Use __builtin_canonicalizef(float) instead of custom flush-to-zero function
Signed-off-by: Aaron Watry <awatry at gmail.com>
CC: Tom Stellard <thomas.stellard at amd.com>
CC: Matt Arsenault <Matthew.Arsenault at amd.com>
---
generic/include/clc/clc.h | 1 +
generic/include/clc/math/remainder.h | 2 +
generic/include/clc/math/remainder.inc | 1 +
generic/lib/SOURCES | 1 +
generic/lib/math/remainder.cl | 508 +++++++++++++++++++++++++++++++++
5 files changed, 513 insertions(+)
create mode 100644 generic/include/clc/math/remainder.h
create mode 100644 generic/include/clc/math/remainder.inc
create mode 100644 generic/lib/math/remainder.cl
diff --git a/generic/include/clc/clc.h b/generic/include/clc/clc.h
index 024bf27..2d4af4b 100644
--- a/generic/include/clc/clc.h
+++ b/generic/include/clc/clc.h
@@ -82,6 +82,7 @@
#include <clc/math/nextafter.h>
#include <clc/math/pow.h>
#include <clc/math/pown.h>
+#include <clc/math/remainder.h>
#include <clc/math/rint.h>
#include <clc/math/round.h>
#include <clc/math/sin.h>
diff --git a/generic/include/clc/math/remainder.h b/generic/include/clc/math/remainder.h
new file mode 100644
index 0000000..97d9fad
--- /dev/null
+++ b/generic/include/clc/math/remainder.h
@@ -0,0 +1,2 @@
+#define __CLC_BODY <clc/math/remainder.inc>
+#include <clc/math/gentype.inc>
diff --git a/generic/include/clc/math/remainder.inc b/generic/include/clc/math/remainder.inc
new file mode 100644
index 0000000..00d0d69
--- /dev/null
+++ b/generic/include/clc/math/remainder.inc
@@ -0,0 +1 @@
+_CLC_OVERLOAD _CLC_DECL __CLC_GENTYPE remainder(__CLC_GENTYPE x, __CLC_GENTYPE y);
diff --git a/generic/lib/SOURCES b/generic/lib/SOURCES
index 517daba..39208bf 100644
--- a/generic/lib/SOURCES
+++ b/generic/lib/SOURCES
@@ -113,6 +113,7 @@ math/tables.cl
math/clc_nextafter.cl
math/nextafter.cl
math/pown.cl
+math/remainder.cl
math/sin.cl
math/sincos.cl
math/sincos_helpers.cl
diff --git a/generic/lib/math/remainder.cl b/generic/lib/math/remainder.cl
new file mode 100644
index 0000000..df42db0
--- /dev/null
+++ b/generic/lib/math/remainder.cl
@@ -0,0 +1,508 @@
+#include <clc/clc.h>
+
+#include "math.h"
+#include "config.h"
+#include "../clcmacro.h"
+
+inline float _clc_remainder_scaleFullRangef32(float y, float t) {
+ float ay, ty, r = 0;
+ int k, iiy, iy, exp_iy0, exp_iy, manty, signy, miy;
+ int delta, shift, ir;
+
+ ay = fabs(t);
+ k = ay > 1024 ? 1024 : (int) ay;
+ k = t < 0 ? -k : k;
+ t = (float) k;
+
+ iiy = as_int(y);
+ iy = iiy & EXSIGNBIT_SP32;
+ signy = iiy & SIGNBIT_SP32;
+ ay = as_float(iy);
+
+ exp_iy0 = iy & EXPBITS_SP32;
+ manty = iy & MANTBITS_SP32;
+
+ //sub-normal
+ ty = exp_iy0 == 0 ? (float) manty : as_float(iy);
+ k = exp_iy0 == 0 ? k - 149 : k;
+ ay = ty;
+ iy = as_int(ay);
+ exp_iy0 = iy & EXPBITS_SP32;
+ exp_iy = (exp_iy0 >> EXPSHIFTBITS_SP32) - EXPBIAS_SP32;
+ // add k to y's exponent
+ r = as_float(iy + (k << EXPSHIFTBITS_SP32));
+ r = (exp_iy + k) > 127 ? as_float(PINFBITPATT_SP32) : r;
+ // add k to y's exponent
+ delta = -126 - (exp_iy + k);
+
+ // sub-normal
+ miy = iy & MANTBITS_SP32;
+ miy |= IMPBIT_SP32;
+ shift = delta > 23 ? 24 : delta;
+ shift = delta < 0 ? 0 : shift;
+ miy >>= shift;
+ r = delta > 0 ? as_float(miy) : r;
+ r = t > (float) (2 * EMAX_SP32) ? as_float(PINFBITPATT_SP32) : r;
+ ir = as_int(r);
+ r = ir <= PINFBITPATT_SP32 ? as_float(as_int(r) | signy) : r;
+ return r;
+}
+/* Scales the float x by 2.0**n.
+Assumes 2*EMIN <= n <= 2*EMAX, though this condition is not checked. */
+inline float _clc_remainder_scaleFloat_2(float x, int n) {
+ float t1, t2;
+ int n1, n2;
+ n1 = n / 2;
+ n2 = n - n1;
+ /* Construct the numbers t1 = 2.0**n1 and t2 = 2.0**n2 */
+ t1 = as_float((n1 + EXPBIAS_SP32) << EXPSHIFTBITS_SP32);
+ t2 = as_float((n2 + EXPBIAS_SP32) << EXPSHIFTBITS_SP32);
+ return (x * t1) * t2;
+}
+/* Scales the float x by 2.0**n.
+ Assumes EMIN <= n <= EMAX, though this condition is not checked. */
+inline float _clc_remainder_scaleFloat_1(float x, int n) {
+ float t;
+ /* Construct the number t = 2.0**n */
+ t = as_float((n + EXPBIAS_SP32) << EXPSHIFTBITS_SP32);
+ return x * t;
+}
+/* Computes the exact product of x and y, the result being the
+nearly double length number (z,zz) */
+inline void _clc_remainder_mul12f(float x, float y, float *z, float *zz) {
+ float hx, tx, hy, ty;
+ // Split x into hx (head) and tx (tail). Do the same for y.
+ uint u;
+ u = as_uint(x);
+ u &= 0xfffff000;
+ hx = as_float(u);
+ tx = x - hx;
+ u = as_uint(y);
+ u &= 0xfffff000;
+ hy = as_float(u);
+ ty = y - hy;
+ *z = x * y;
+ *zz = (((hx * hy - *z) + hx * ty) + tx * hy) + tx * ty;
+}
+
+_CLC_OVERLOAD _CLC_DEF float remainder(float x, float y) {
+ if (!__clc_fp32_subnormals_supported()) {
+ const int loop_scale = 12;
+ const float fscale = 1.0f / (float) (1 << loop_scale);
+
+ int ntimes;
+ float ret = 0;
+ int ui_x, ui_y, ui_ax, ui_ay, xexp, yexp, signx;
+ float af_x, af_y, af_ybase, fx, fxp, fxm, fy, w, scale, t, c, cc, v;
+ float yscale, scaled_w, saved_w, div, sdiv, ratio, sratio, fxexp, sub_fx;
+ int iw_scaled, wexp, it, i, ifx, ex, ey;
+ ;
+ float xr, xr0, xr_base, yr;
+ uint q;
+
+ ui_x = as_int(x);
+ ui_y = as_int(y);
+ ui_ax = ui_x & EXSIGNBIT_SP32;
+ ui_ay = ui_y & EXSIGNBIT_SP32;
+
+ /* special case handle */
+ if (ui_ax > PINFBITPATT_SP32)
+ return x;
+ if (ui_ax == PINFBITPATT_SP32)
+ return as_float(QNANBITPATT_SP32);
+ if (ui_ay > PINFBITPATT_SP32)
+ return y;
+ if (ui_ay == PINFBITPATT_SP32)
+ return x;
+ if (ui_ay == 0 && ui_ax == 0)
+ return as_float(QNANBITPATT_SP32);
+ if (ui_ax == 0)
+ return x;
+ if (ui_ay == 0)
+ return as_float(QNANBITPATT_SP32);
+
+ signx = ui_x & SIGNBIT_SP32;
+ af_x = as_float(ui_ax);
+ af_ybase = af_y = as_float(ui_ay);
+ yexp = (int) ((ui_y & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);
+
+ yscale = (float) ((yexp < 48 && ui_ay != 0) ? (48 - yexp) : 0);
+ if (yscale != 0) {
+ af_y = _clc_remainder_scaleFullRangef32(af_ybase, yscale);
+ }
+
+ ui_y = as_int(af_y);
+ yexp = (int) ((ui_y & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);
+ xexp = (int) ((ui_x & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);
+ fx = af_x;
+ fy = af_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 24 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. */
+ ntimes = (xexp - yexp) / loop_scale;
+ ntimes = xexp <= yexp ? 0 : ntimes;
+
+ /* Set w = y * 2^(ntimes*loop_scale) */
+ w = _clc_remainder_scaleFloat_2(fy, ntimes * loop_scale);
+ w = ntimes == 0 ? fy : w;
+
+ /* Set scale = 2^(-loop_scale) */
+ scale = ntimes == 0 ? 1.0f : fscale;
+
+ // make sure recip does not overflow
+ wexp = (int) ((as_int(w) & EXPBITS_SP32) >> EXPSHIFTBITS_SP32) - EXPBIAS_SP32;
+ saved_w = w;
+ scaled_w = _clc_remainder_scaleFloat_1(w, -14);
+ iw_scaled = wexp > 105 & wexp <= 127;
+ w = (iw_scaled & ntimes) > 0 ? scaled_w : w;
+
+ /* 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 can be performed exactly when performed
+ in double precision, and the result at each stage can
+ fit exactly in a single precision number. */
+ for (i = 0; i < ntimes; i++) {
+ /* Set fx = fx - w * t, where t is equal to trunc(dx/w). */
+ div = fx / w; //was __amdil_improved_div_f32 => native_div(x, y) => native_recip(y)*x
+ sdiv = _clc_remainder_scaleFloat_1(div, -14);
+ div = iw_scaled ? sdiv : div;
+ t = floor(div);
+ w = saved_w;
+ iw_scaled = 0;
+
+ /* At this point, t may be one too large due to rounding of fx/w */
+
+ /* Compute w * t in quad precision */
+ _clc_remainder_mul12f(w, t, &c, &cc);
+
+ /* Subtract w * t from fx */
+ v = fx - c;
+ fx = v + (((fx - v) - c) - cc);
+
+ /* If t was one too large, fx will be negative. Add back one w */
+ /* It might be possible to speed up this loop by finding
+ a way to compute correctly truncated t directly from fx and w.
+ We would then avoid the need for this check on negative fx. */
+ fxp = fx + w;
+ fxm = fx - w;
+ fx = fx < 0.0f ? fxp : fx;
+ fx = fx >= w ? fxm : fx;
+
+ /* Scale w down by for the next iteration */
+ w *= scale;
+ saved_w = w;
+ }
+
+ /* One more time */
+ // iw = as_int(w);
+ ifx = as_int(fx);
+ fxexp = (int) ((ifx & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);
+ // wexp = (int) ((iw & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);
+ sub_fx = fx;
+ // make sure recip does not overflow
+ wexp = (int) ((as_int(w) & EXPBITS_SP32) >> EXPSHIFTBITS_SP32) - EXPBIAS_SP32;
+ saved_w = w;
+ scaled_w = _clc_remainder_scaleFloat_1(w, -14);
+ iw_scaled = wexp > 105 & wexp <= 127;
+ w = iw_scaled ? scaled_w : w;
+ ratio = fx / w; //was the amdil equivalent of native_divide(x, y) => native_recip(y)*x;
+ sratio = _clc_remainder_scaleFloat_1(ratio, -14);
+ ratio = iw_scaled ? sratio : ratio;
+ t = floor(ratio);
+ it = (int) t;
+
+ w = saved_w;
+ _clc_remainder_mul12f(w, t, &c, &cc);
+
+ v = fx - c;
+ fx = v + (((fx - v) - c) - cc);
+
+ if (fx < 0.0f) {
+ fx += w;
+ it--;
+ }
+
+ if (fx >= w) {
+ fx -= w;
+ it++;
+ }
+
+ // sub-normal fax
+ fx = fxexp == 0 ? sub_fx : fx;
+
+ float scaleback = 0;
+
+ // in case fx == 0 and we'got a divisor
+ it = (yscale > 30) ? 0 : ((unsigned int) it << (int) yscale);
+
+ if (as_int(fx) != 0 && yscale != 0) {
+ xr = fx;
+ xr_base = fx;
+ yr = af_ybase;
+ q = 0;
+ ex = ilogb(fx);
+ ey = ilogb(af_ybase);
+
+ yr = (float) _clc_remainder_scaleFullRangef32(af_ybase, (float) -ey);
+ xr = (float) _clc_remainder_scaleFullRangef32(fx, (float) -ex);
+
+ for (i = ex - ey; i > 0; i--) {
+ q <<= 1;
+ xr0 = xr;
+ xr = (xr0 >= yr) ? xr0 - yr : xr0;
+ q = (xr0 >= yr) ? q + 1 : q;
+ xr += xr;
+ }
+ q <<= 1;
+ xr0 = xr;
+ xr = (xr0 >= yr) ? xr0 - yr : xr0;
+ q = (xr0 >= yr) ? q + 1 : q;
+ xr = _clc_remainder_scaleFullRangef32(xr, (float) ey);
+
+ fx = (ex - ey >= 0) ? xr : xr_base;
+ q = (ex - ey >= 0) ? q : 0;
+ it += q;
+
+ xexp = (int) ((as_int(fx) & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);
+
+ w = af_ybase;
+ if (xexp < 24) {
+ fx = _clc_remainder_scaleFullRangef32(fx, 48);
+ w = _clc_remainder_scaleFullRangef32(af_ybase, 48);
+ scaleback = -48;
+ }
+ }
+ /* 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 (== fy == y) if necessary. */
+ if (fx * 2.f > w || ((fx * 2.f == w) && (it & 1))) {
+ fx -= w;
+ it++;
+ }
+ if (scaleback != 0) {
+ fx = _clc_remainder_scaleFullRangef32(fx, scaleback);
+ }
+
+ ret = (signx) ? as_float(as_int(fx) ^ SIGNBIT_SP32) : fx;
+
+ return ret;
+
+ }
+
+ //Otherwise, Subnormals are supported
+
+ x = __builtin_canonicalizef(x);
+ y = __builtin_canonicalizef(y);
+
+ 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 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;
+
+#define _CLC_BIT c = xr >= yr; q = (q << 1) | c; xr -= c ? yr : 0.0f; xr += xr
+
+ while (k > 3) {
+ _CLC_BIT;
+ _CLC_BIT;
+ _CLC_BIT;
+ _CLC_BIT;
+ k -= 4;
+ }
+
+ while (k > 0) {
+ _CLC_BIT;
+ --k;
+ }
+
+#undef _CLC_BIT
+
+ 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;
+
+ c = ax == ay;
+ xr = c ? 0.0f : xr;
+
+ xr = as_float(sx ^ as_int(xr));
+
+ c = ax > PINFBITPATT_SP32 | ay > PINFBITPATT_SP32 | ax == PINFBITPATT_SP32 | ay == 0;
+ xr = c ? as_float(QNANBITPATT_SP32) : xr;
+
+ return xr;
+}
+
+_CLC_BINARY_VECTORIZE(_CLC_OVERLOAD _CLC_DEF, float, remainder, float, float);
+
+#ifdef cl_khr_fp64
+
+#pragma OPENCL EXTENSION cl_khr_fp64 : enable
+
+inline double
+_clc_remainder_ldexp(double x, int n) {
+ // XXX Have to go twice here because the hardware can't handle the full range (yet)
+ int nh = n >> 1;
+ return ldexp(ldexp(x, nh), n - nh);
+}
+
+_CLC_OVERLOAD _CLC_DEF double remainder(double y, double x) {
+ 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 = _clc_remainder_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;
+
+ // At this point, dx lies in the range [0,dy)
+
+ // For the fmod function, we're done apart from setting the correct sign.
+ //
+ // 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;
+
+ double ret = as_double(xsgn ^ as_ulong(dx));
+ dx = as_double(ax);
+
+ // Now handle |x| == |y|
+ int c = dx == dy;
+ t = as_double(xsgn);
+ ret = c ? t : ret;
+
+ // Next, handle |x| < |y|
+ c = dx < dy;
+ ret = c ? x : ret;
+
+ c &= (yexp < 1023 & 2.0 * dx > dy) | (dx > 0.5 * dy);
+ // 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;
+ ret = c ? as_double(QNANBITPATT_DP64) : ret;
+
+ // y is +-Inf, NaN
+ c = yexp > BIASEDEMAX_DP64;
+ t = y == y ? x : y;
+ ret = c ? t : ret;
+
+ // x is +=Inf, NaN
+ c = xexp > BIASEDEMAX_DP64;
+ ret = c ? as_double(QNANBITPATT_DP64) : ret;
+
+ return ret;
+}
+
+_CLC_BINARY_VECTORIZE(_CLC_OVERLOAD _CLC_DEF, double, remainder, double, double);
+#endif
\ No newline at end of file
--
2.9.3
More information about the Libclc-dev
mailing list