<html><head><meta http-equiv="Content-Type" content="text/html charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space;" class=""><br class=""><div><blockquote type="cite" class=""><div class="">On Feb 6, 2017, at 20:32, Aaron Watry <<a href="mailto:awatry@gmail.com" class="">awatry@gmail.com</a>> wrote:</div><br class="Apple-interchange-newline"><div class=""><div dir="ltr" class=""><br class=""><div class="gmail_extra"><br class=""><div class="gmail_quote">On Thu, Jan 19, 2017 at 2:47 PM, Matt Arsenault <span dir="ltr" class=""><<a href="mailto:arsenm2@gmail.com" target="_blank" class="">arsenm2@gmail.com</a>></span> wrote:<br class=""><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">This fails conformance for me:<br class="">
<br class="">
<br class=""></blockquote><div class=""><br class="">Interestingly enough some (but not all) of the test inputs start to pass when I forcefully enable subnormal support in libclc (--enable-runtime-subnormal and also updating generic/lib/shared/<a href="http://subnormal_config.cl/" class="">subnormal_config.cl</a> to enable 32-bit subnormals).<br class=""><br class=""></div><div class="">I'll be playing with this a bit as time permits.<br class=""><br class=""></div><div class="">—Aaron<br class=""></div></div></div></div></div></blockquote><div><br class=""></div><div>I’m guessing you need to somehow access one of the other division implementations. By default fdiv will be getting the !fpmath 2.5 ULP metadata. Can you try adding -cl-f432-correctly-rounded-divide-sqrt to the build of this file, or calling a wrapper IR function which avoids the metadata?</div><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div class="gmail_extra"><div class="gmail_quote"><div class=""><br class=""> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br class="">
> On Jan 18, 2017, at 20:04, Aaron Watry via Libclc-dev <<a href="mailto:libclc-dev@lists.llvm.org" class="">libclc-dev@lists.llvm.org</a>> wrote:<br class="">
><br class="">
> Mostly ported from the amd-builtins branch.<br class="">
><br class="">
> The amd-builtins branch uses __amdil_improved_fdiv_f32 and FTZ which aren't available in generic CLC.<br class="">
><br class="">
> __amdil_improved_fdiv_f32 points to native_divide which does native_recip(y)*x.<br class="">
><br class="">
> Since we don't have native_divide or native_recip yet, I've just stuck an actual division here.<br class="">
><br class="">
> I've taken a shot at a replacement for FTZ(x), but feel free to suggest alternatives.<br class="">
><br class="">
> Tested via piglit on a Radeon HD 7850 using the tests just sent to that list.<br class="">
><br class="">
> v2: Use __builtin_canonicalizef(float) instead of custom flush-to-zero function<br class="">
><br class="">
> Signed-off-by: Aaron Watry <<a href="mailto:awatry@gmail.com" class="">awatry@gmail.com</a>><br class="">
> CC: Tom Stellard <<a href="mailto:thomas.stellard@amd.com" class="">thomas.stellard@amd.com</a>><br class="">
> CC: Matt Arsenault <<a href="mailto:Matthew.Arsenault@amd.com" class="">Matthew.Arsenault@amd.com</a>><br class="">
> ---<br class="">
> generic/include/clc/clc.h | 1 +<br class="">
> generic/include/clc/math/<wbr class="">remainder.h | 2 +<br class="">
> generic/include/clc/math/<wbr class="">remainder.inc | 1 +<br class="">
> generic/lib/SOURCES | 1 +<br class="">
> generic/lib/math/<a href="http://remainder.cl/" rel="noreferrer" target="_blank" class="">remainder.cl</a> | 508 ++++++++++++++++++++++++++++++<wbr class="">+++<br class="">
> 5 files changed, 513 insertions(+)<br class="">
> create mode 100644 generic/include/clc/math/<wbr class="">remainder.h<br class="">
> create mode 100644 generic/include/clc/math/<wbr class="">remainder.inc<br class="">
> create mode 100644 generic/lib/math/<a href="http://remainder.cl/" rel="noreferrer" target="_blank" class="">remainder.cl</a><br class="">
><br class="">
> diff --git a/generic/include/clc/clc.h b/generic/include/clc/clc.h<br class="">
> index 024bf27..2d4af4b 100644<br class="">
> --- a/generic/include/clc/clc.h<br class="">
> +++ b/generic/include/clc/clc.h<br class="">
> @@ -82,6 +82,7 @@<br class="">
> #include <clc/math/nextafter.h><br class="">
> #include <clc/math/pow.h><br class="">
> #include <clc/math/pown.h><br class="">
> +#include <clc/math/remainder.h><br class="">
> #include <clc/math/rint.h><br class="">
> #include <clc/math/round.h><br class="">
> #include <clc/math/sin.h><br class="">
> diff --git a/generic/include/clc/math/<wbr class="">remainder.h b/generic/include/clc/math/<wbr class="">remainder.h<br class="">
> new file mode 100644<br class="">
> index 0000000..97d9fad<br class="">
> --- /dev/null<br class="">
> +++ b/generic/include/clc/math/<wbr class="">remainder.h<br class="">
> @@ -0,0 +1,2 @@<br class="">
> +#define __CLC_BODY <clc/math/remainder.inc><br class="">
> +#include <clc/math/gentype.inc><br class="">
> diff --git a/generic/include/clc/math/<wbr class="">remainder.inc b/generic/include/clc/math/<wbr class="">remainder.inc<br class="">
> new file mode 100644<br class="">
> index 0000000..00d0d69<br class="">
> --- /dev/null<br class="">
> +++ b/generic/include/clc/math/<wbr class="">remainder.inc<br class="">
> @@ -0,0 +1 @@<br class="">
> +_CLC_OVERLOAD _CLC_DECL __CLC_GENTYPE remainder(__CLC_GENTYPE x, __CLC_GENTYPE y);<br class="">
> diff --git a/generic/lib/SOURCES b/generic/lib/SOURCES<br class="">
> index 517daba..39208bf 100644<br class="">
> --- a/generic/lib/SOURCES<br class="">
> +++ b/generic/lib/SOURCES<br class="">
> @@ -113,6 +113,7 @@ math/<a href="http://tables.cl/" rel="noreferrer" target="_blank" class="">tables.cl</a><br class="">
> math/<a href="http://clc_nextafter.cl/" rel="noreferrer" target="_blank" class="">clc_nextafter.cl</a><br class="">
> math/<a href="http://nextafter.cl/" rel="noreferrer" target="_blank" class="">nextafter.cl</a><br class="">
> math/<a href="http://pown.cl/" rel="noreferrer" target="_blank" class="">pown.cl</a><br class="">
> +math/<a href="http://remainder.cl/" rel="noreferrer" target="_blank" class="">remainder.cl</a><br class="">
> math/<a href="http://sin.cl/" rel="noreferrer" target="_blank" class="">sin.cl</a><br class="">
> math/<a href="http://sincos.cl/" rel="noreferrer" target="_blank" class="">sincos.cl</a><br class="">
> math/<a href="http://sincos_helpers.cl/" rel="noreferrer" target="_blank" class="">sincos_helpers.cl</a><br class="">
> diff --git a/generic/lib/math/<a href="http://remainder.cl/" rel="noreferrer" target="_blank" class="">remainder.<wbr class="">cl</a> b/generic/lib/math/<a href="http://remainder.cl/" rel="noreferrer" target="_blank" class="">remainder.<wbr class="">cl</a><br class="">
> new file mode 100644<br class="">
> index 0000000..df42db0<br class="">
> --- /dev/null<br class="">
> +++ b/generic/lib/math/<a href="http://remainder.cl/" rel="noreferrer" target="_blank" class="">remainder.<wbr class="">cl</a><br class="">
> @@ -0,0 +1,508 @@<br class="">
> +#include <clc/clc.h><br class="">
> +<br class="">
> +#include "math.h"<br class="">
> +#include "config.h"<br class="">
> +#include "../clcmacro.h"<br class="">
> +<br class="">
> +inline float _clc_remainder_<wbr class="">scaleFullRangef32(float y, float t) {<br class="">
> + float ay, ty, r = 0;<br class="">
> + int k, iiy, iy, exp_iy0, exp_iy, manty, signy, miy;<br class="">
> + int delta, shift, ir;<br class="">
> +<br class="">
> + ay = fabs(t);<br class="">
> + k = ay > 1024 ? 1024 : (int) ay;<br class="">
> + k = t < 0 ? -k : k;<br class="">
> + t = (float) k;<br class="">
> +<br class="">
> + iiy = as_int(y);<br class="">
> + iy = iiy & EXSIGNBIT_SP32;<br class="">
> + signy = iiy & SIGNBIT_SP32;<br class="">
> + ay = as_float(iy);<br class="">
> +<br class="">
> + exp_iy0 = iy & EXPBITS_SP32;<br class="">
> + manty = iy & MANTBITS_SP32;<br class="">
> +<br class="">
> + //sub-normal<br class="">
> + ty = exp_iy0 == 0 ? (float) manty : as_float(iy);<br class="">
> + k = exp_iy0 == 0 ? k - 149 : k;<br class="">
> + ay = ty;<br class="">
> + iy = as_int(ay);<br class="">
> + exp_iy0 = iy & EXPBITS_SP32;<br class="">
> + exp_iy = (exp_iy0 >> EXPSHIFTBITS_SP32) - EXPBIAS_SP32;<br class="">
> + // add k to y's exponent<br class="">
> + r = as_float(iy + (k << EXPSHIFTBITS_SP32));<br class="">
> + r = (exp_iy + k) > 127 ? as_float(PINFBITPATT_SP32) : r;<br class="">
> + // add k to y's exponent<br class="">
> + delta = -126 - (exp_iy + k);<br class="">
> +<br class="">
> + // sub-normal<br class="">
> + miy = iy & MANTBITS_SP32;<br class="">
> + miy |= IMPBIT_SP32;<br class="">
> + shift = delta > 23 ? 24 : delta;<br class="">
> + shift = delta < 0 ? 0 : shift;<br class="">
> + miy >>= shift;<br class="">
> + r = delta > 0 ? as_float(miy) : r;<br class="">
> + r = t > (float) (2 * EMAX_SP32) ? as_float(PINFBITPATT_SP32) : r;<br class="">
> + ir = as_int(r);<br class="">
> + r = ir <= PINFBITPATT_SP32 ? as_float(as_int(r) | signy) : r;<br class="">
> + return r;<br class="">
> +}<br class="">
> +/* Scales the float x by 2.0**n.<br class="">
> +Assumes 2*EMIN <= n <= 2*EMAX, though this condition is not checked. */<br class="">
> +inline float _clc_remainder_scaleFloat_2(<wbr class="">float x, int n) {<br class="">
> + float t1, t2;<br class="">
> + int n1, n2;<br class="">
> + n1 = n / 2;<br class="">
> + n2 = n - n1;<br class="">
> + /* Construct the numbers t1 = 2.0**n1 and t2 = 2.0**n2 */<br class="">
> + t1 = as_float((n1 + EXPBIAS_SP32) << EXPSHIFTBITS_SP32);<br class="">
> + t2 = as_float((n2 + EXPBIAS_SP32) << EXPSHIFTBITS_SP32);<br class="">
> + return (x * t1) * t2;<br class="">
> +}<br class="">
> +/* Scales the float x by 2.0**n.<br class="">
> + Assumes EMIN <= n <= EMAX, though this condition is not checked. */<br class="">
> +inline float _clc_remainder_scaleFloat_1(<wbr class="">float x, int n) {<br class="">
> + float t;<br class="">
> + /* Construct the number t = 2.0**n */<br class="">
> + t = as_float((n + EXPBIAS_SP32) << EXPSHIFTBITS_SP32);<br class="">
> + return x * t;<br class="">
> +}<br class="">
> +/* Computes the exact product of x and y, the result being the<br class="">
> +nearly double length number (z,zz) */<br class="">
> +inline void _clc_remainder_mul12f(float x, float y, float *z, float *zz) {<br class="">
> + float hx, tx, hy, ty;<br class="">
> + // Split x into hx (head) and tx (tail). Do the same for y.<br class="">
> + uint u;<br class="">
> + u = as_uint(x);<br class="">
> + u &= 0xfffff000;<br class="">
> + hx = as_float(u);<br class="">
> + tx = x - hx;<br class="">
> + u = as_uint(y);<br class="">
> + u &= 0xfffff000;<br class="">
> + hy = as_float(u);<br class="">
> + ty = y - hy;<br class="">
> + *z = x * y;<br class="">
> + *zz = (((hx * hy - *z) + hx * ty) + tx * hy) + tx * ty;<br class="">
> +}<br class="">
> +<br class="">
> +_CLC_OVERLOAD _CLC_DEF float remainder(float x, float y) {<br class="">
> + if (!__clc_fp32_subnormals_<wbr class="">supported()) {<br class="">
> + const int loop_scale = 12;<br class="">
> + const float fscale = 1.0f / (float) (1 << loop_scale);<br class="">
> +<br class="">
> + int ntimes;<br class="">
> + float ret = 0;<br class="">
> + int ui_x, ui_y, ui_ax, ui_ay, xexp, yexp, signx;<br class="">
> + float af_x, af_y, af_ybase, fx, fxp, fxm, fy, w, scale, t, c, cc, v;<br class="">
> + float yscale, scaled_w, saved_w, div, sdiv, ratio, sratio, fxexp, sub_fx;<br class="">
> + int iw_scaled, wexp, it, i, ifx, ex, ey;<br class="">
> + ;<br class="">
> + float xr, xr0, xr_base, yr;<br class="">
> + uint q;<br class="">
> +<br class="">
> + ui_x = as_int(x);<br class="">
> + ui_y = as_int(y);<br class="">
> + ui_ax = ui_x & EXSIGNBIT_SP32;<br class="">
> + ui_ay = ui_y & EXSIGNBIT_SP32;<br class="">
> +<br class="">
> + /* special case handle */<br class="">
> + if (ui_ax > PINFBITPATT_SP32)<br class="">
> + return x;<br class="">
> + if (ui_ax == PINFBITPATT_SP32)<br class="">
> + return as_float(QNANBITPATT_SP32);<br class="">
> + if (ui_ay > PINFBITPATT_SP32)<br class="">
> + return y;<br class="">
> + if (ui_ay == PINFBITPATT_SP32)<br class="">
> + return x;<br class="">
> + if (ui_ay == 0 && ui_ax == 0)<br class="">
> + return as_float(QNANBITPATT_SP32);<br class="">
> + if (ui_ax == 0)<br class="">
> + return x;<br class="">
> + if (ui_ay == 0)<br class="">
> + return as_float(QNANBITPATT_SP32);<br class="">
> +<br class="">
> + signx = ui_x & SIGNBIT_SP32;<br class="">
> + af_x = as_float(ui_ax);<br class="">
> + af_ybase = af_y = as_float(ui_ay);<br class="">
> + yexp = (int) ((ui_y & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);<br class="">
> +<br class="">
> + yscale = (float) ((yexp < 48 && ui_ay != 0) ? (48 - yexp) : 0);<br class="">
> + if (yscale != 0) {<br class="">
> + af_y = _clc_remainder_<wbr class="">scaleFullRangef32(af_ybase, yscale);<br class="">
> + }<br class="">
> +<br class="">
> + ui_y = as_int(af_y);<br class="">
> + yexp = (int) ((ui_y & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);<br class="">
> + xexp = (int) ((ui_x & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);<br class="">
> + fx = af_x;<br class="">
> + fy = af_y;<br class="">
> +<br class="">
> + /* Set ntimes to the number of times we need to do a<br class="">
> + partial remainder. If the exponent of x is an exact multiple<br class="">
> + of 24 larger than the exponent of y, and the mantissa of x is<br class="">
> + less than the mantissa of y, ntimes will be one too large<br class="">
> + but it doesn't matter - it just means that we'll go round<br class="">
> + the loop below one extra time. */<br class="">
> + ntimes = (xexp - yexp) / loop_scale;<br class="">
> + ntimes = xexp <= yexp ? 0 : ntimes;<br class="">
> +<br class="">
> + /* Set w = y * 2^(ntimes*loop_scale) */<br class="">
> + w = _clc_remainder_scaleFloat_2(<wbr class="">fy, ntimes * loop_scale);<br class="">
> + w = ntimes == 0 ? fy : w;<br class="">
> +<br class="">
> + /* Set scale = 2^(-loop_scale) */<br class="">
> + scale = ntimes == 0 ? 1.0f : fscale;<br class="">
> +<br class="">
> + // make sure recip does not overflow<br class="">
> + wexp = (int) ((as_int(w) & EXPBITS_SP32) >> EXPSHIFTBITS_SP32) - EXPBIAS_SP32;<br class="">
> + saved_w = w;<br class="">
> + scaled_w = _clc_remainder_scaleFloat_1(w, -14);<br class="">
> + iw_scaled = wexp > 105 & wexp <= 127;<br class="">
> + w = (iw_scaled & ntimes) > 0 ? scaled_w : w;<br class="">
> +<br class="">
> + /* Each time round the loop we compute a partial remainder.<br class="">
> + This is done by subtracting a large multiple of w<br class="">
> + from x each time, where w is a scaled up version of y.<br class="">
> + The subtraction can be performed exactly when performed<br class="">
> + in double precision, and the result at each stage can<br class="">
> + fit exactly in a single precision number. */<br class="">
> + for (i = 0; i < ntimes; i++) {<br class="">
> + /* Set fx = fx - w * t, where t is equal to trunc(dx/w). */<br class="">
> + div = fx / w; //was __amdil_improved_div_f32 => native_div(x, y) => native_recip(y)*x<br class="">
> + sdiv = _clc_remainder_scaleFloat_1(<wbr class="">div, -14);<br class="">
> + div = iw_scaled ? sdiv : div;<br class="">
> + t = floor(div);<br class="">
> + w = saved_w;<br class="">
> + iw_scaled = 0;<br class="">
> +<br class="">
> + /* At this point, t may be one too large due to rounding of fx/w */<br class="">
> +<br class="">
> + /* Compute w * t in quad precision */<br class="">
> + _clc_remainder_mul12f(w, t, &c, &cc);<br class="">
> +<br class="">
> + /* Subtract w * t from fx */<br class="">
> + v = fx - c;<br class="">
> + fx = v + (((fx - v) - c) - cc);<br class="">
> +<br class="">
> + /* If t was one too large, fx will be negative. Add back one w */<br class="">
> + /* It might be possible to speed up this loop by finding<br class="">
> + a way to compute correctly truncated t directly from fx and w.<br class="">
> + We would then avoid the need for this check on negative fx. */<br class="">
> + fxp = fx + w;<br class="">
> + fxm = fx - w;<br class="">
> + fx = fx < 0.0f ? fxp : fx;<br class="">
> + fx = fx >= w ? fxm : fx;<br class="">
> +<br class="">
> + /* Scale w down by for the next iteration */<br class="">
> + w *= scale;<br class="">
> + saved_w = w;<br class="">
> + }<br class="">
> +<br class="">
> + /* One more time */<br class="">
> + // iw = as_int(w);<br class="">
> + ifx = as_int(fx);<br class="">
> + fxexp = (int) ((ifx & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);<br class="">
> + // wexp = (int) ((iw & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);<br class="">
> + sub_fx = fx;<br class="">
> + // make sure recip does not overflow<br class="">
> + wexp = (int) ((as_int(w) & EXPBITS_SP32) >> EXPSHIFTBITS_SP32) - EXPBIAS_SP32;<br class="">
> + saved_w = w;<br class="">
> + scaled_w = _clc_remainder_scaleFloat_1(w, -14);<br class="">
> + iw_scaled = wexp > 105 & wexp <= 127;<br class="">
> + w = iw_scaled ? scaled_w : w;<br class="">
> + ratio = fx / w; //was the amdil equivalent of native_divide(x, y) => native_recip(y)*x;<br class="">
> + sratio = _clc_remainder_scaleFloat_1(<wbr class="">ratio, -14);<br class="">
> + ratio = iw_scaled ? sratio : ratio;<br class="">
> + t = floor(ratio);<br class="">
> + it = (int) t;<br class="">
> +<br class="">
> + w = saved_w;<br class="">
> + _clc_remainder_mul12f(w, t, &c, &cc);<br class="">
> +<br class="">
> + v = fx - c;<br class="">
> + fx = v + (((fx - v) - c) - cc);<br class="">
> +<br class="">
> + if (fx < 0.0f) {<br class="">
> + fx += w;<br class="">
> + it--;<br class="">
> + }<br class="">
> +<br class="">
> + if (fx >= w) {<br class="">
> + fx -= w;<br class="">
> + it++;<br class="">
> + }<br class="">
> +<br class="">
> + // sub-normal fax<br class="">
> + fx = fxexp == 0 ? sub_fx : fx;<br class="">
> +<br class="">
> + float scaleback = 0;<br class="">
> +<br class="">
> + // in case fx == 0 and we'got a divisor<br class="">
> + it = (yscale > 30) ? 0 : ((unsigned int) it << (int) yscale);<br class="">
> +<br class="">
> + if (as_int(fx) != 0 && yscale != 0) {<br class="">
> + xr = fx;<br class="">
> + xr_base = fx;<br class="">
> + yr = af_ybase;<br class="">
> + q = 0;<br class="">
> + ex = ilogb(fx);<br class="">
> + ey = ilogb(af_ybase);<br class="">
> +<br class="">
> + yr = (float) _clc_remainder_<wbr class="">scaleFullRangef32(af_ybase, (float) -ey);<br class="">
> + xr = (float) _clc_remainder_<wbr class="">scaleFullRangef32(fx, (float) -ex);<br class="">
> +<br class="">
> + for (i = ex - ey; i > 0; i--) {<br class="">
> + q <<= 1;<br class="">
> + xr0 = xr;<br class="">
> + xr = (xr0 >= yr) ? xr0 - yr : xr0;<br class="">
> + q = (xr0 >= yr) ? q + 1 : q;<br class="">
> + xr += xr;<br class="">
> + }<br class="">
> + q <<= 1;<br class="">
> + xr0 = xr;<br class="">
> + xr = (xr0 >= yr) ? xr0 - yr : xr0;<br class="">
> + q = (xr0 >= yr) ? q + 1 : q;<br class="">
> + xr = _clc_remainder_<wbr class="">scaleFullRangef32(xr, (float) ey);<br class="">
> +<br class="">
> + fx = (ex - ey >= 0) ? xr : xr_base;<br class="">
> + q = (ex - ey >= 0) ? q : 0;<br class="">
> + it += q;<br class="">
> +<br class="">
> + xexp = (int) ((as_int(fx) & EXPBITS_SP32) >> EXPSHIFTBITS_SP32);<br class="">
> +<br class="">
> + w = af_ybase;<br class="">
> + if (xexp < 24) {<br class="">
> + fx = _clc_remainder_<wbr class="">scaleFullRangef32(fx, 48);<br class="">
> + w = _clc_remainder_<wbr class="">scaleFullRangef32(af_ybase, 48);<br class="">
> + scaleback = -48;<br class="">
> + }<br class="">
> + }<br class="">
> + /* At this point, dx lies in the range [0,dy) */<br class="">
> + /* For the remainder function, we need to adjust dx<br class="">
> + so that it lies in the range (-y/2, y/2] by carefully<br class="">
> + subtracting w (== fy == y) if necessary. */<br class="">
> + if (fx * 2.f > w || ((fx * 2.f == w) && (it & 1))) {<br class="">
> + fx -= w;<br class="">
> + it++;<br class="">
> + }<br class="">
> + if (scaleback != 0) {<br class="">
> + fx = _clc_remainder_<wbr class="">scaleFullRangef32(fx, scaleback);<br class="">
> + }<br class="">
> +<br class="">
> + ret = (signx) ? as_float(as_int(fx) ^ SIGNBIT_SP32) : fx;<br class="">
> +<br class="">
> + return ret;<br class="">
> +<br class="">
> + }<br class="">
> +<br class="">
> + //Otherwise, Subnormals are supported<br class="">
> +<br class="">
> + x = __builtin_canonicalizef(x);<br class="">
> + y = __builtin_canonicalizef(y);<br class="">
> +<br class="">
> + int ux = as_int(x);<br class="">
> + int ax = ux & EXSIGNBIT_SP32;<br class="">
> + float xa = as_float(ax);<br class="">
> + int sx = ux ^ ax;<br class="">
> + int ex = ax >> EXPSHIFTBITS_SP32;<br class="">
> +<br class="">
> + int uy = as_int(y);<br class="">
> + int ay = uy & EXSIGNBIT_SP32;<br class="">
> + float ya = as_float(ay);<br class="">
> + int ey = ay >> EXPSHIFTBITS_SP32;<br class="">
> +<br class="">
> + float xr = as_float(0x3f800000 | (ax & 0x007fffff));<br class="">
> + float yr = as_float(0x3f800000 | (ay & 0x007fffff));<br class="">
> + int c;<br class="">
> + int k = ex - ey;<br class="">
> +<br class="">
> + uint q = 0;<br class="">
> +<br class="">
> +#define _CLC_BIT c = xr >= yr; q = (q << 1) | c; xr -= c ? yr : 0.0f; xr += xr<br class="">
> +<br class="">
> + while (k > 3) {<br class="">
> + _CLC_BIT;<br class="">
> + _CLC_BIT;<br class="">
> + _CLC_BIT;<br class="">
> + _CLC_BIT;<br class="">
> + k -= 4;<br class="">
> + }<br class="">
> +<br class="">
> + while (k > 0) {<br class="">
> + _CLC_BIT;<br class="">
> + --k;<br class="">
> + }<br class="">
> +<br class="">
> +#undef _CLC_BIT<br class="">
> +<br class="">
> + c = xr > yr;<br class="">
> + q = (q << 1) | c;<br class="">
> + xr -= c ? yr : 0.0f;<br class="">
> +<br class="">
> + int lt = ex < ey;<br class="">
> +<br class="">
> + q = lt ? 0 : q;<br class="">
> + xr = lt ? xa : xr;<br class="">
> + yr = lt ? ya : yr;<br class="">
> +<br class="">
> + c = (yr < 2.0f * xr) | ((yr == 2.0f * xr) & ((q & 0x1) == 0x1));<br class="">
> + xr -= c ? yr : 0.0f;<br class="">
> + q += c;<br class="">
> +<br class="">
> + float s = as_float(ey << EXPSHIFTBITS_SP32);<br class="">
> + xr *= lt ? 1.0f : s;<br class="">
> +<br class="">
> + c = ax == ay;<br class="">
> + xr = c ? 0.0f : xr;<br class="">
> +<br class="">
> + xr = as_float(sx ^ as_int(xr));<br class="">
> +<br class="">
> + c = ax > PINFBITPATT_SP32 | ay > PINFBITPATT_SP32 | ax == PINFBITPATT_SP32 | ay == 0;<br class="">
> + xr = c ? as_float(QNANBITPATT_SP32) : xr;<br class="">
> +<br class="">
> + return xr;<br class="">
> +}<br class="">
> +<br class="">
> +_CLC_BINARY_VECTORIZE(_CLC_<wbr class="">OVERLOAD _CLC_DEF, float, remainder, float, float);<br class="">
> +<br class="">
> +#ifdef cl_khr_fp64<br class="">
> +<br class="">
> +#pragma OPENCL EXTENSION cl_khr_fp64 : enable<br class="">
> +<br class="">
> +inline double<br class="">
> +_clc_remainder_ldexp(double x, int n) {<br class="">
> + // XXX Have to go twice here because the hardware can't handle the full range (yet)<br class="">
> + int nh = n >> 1;<br class="">
> + return ldexp(ldexp(x, nh), n - nh);<br class="">
> +}<br class="">
> +<br class="">
> +_CLC_OVERLOAD _CLC_DEF double remainder(double y, double x) {<br class="">
> + ulong ux = as_ulong(x);<br class="">
> + ulong ax = ux & ~SIGNBIT_DP64;<br class="">
> + ulong xsgn = ux ^ ax;<br class="">
> + double dx = as_double(ax);<br class="">
> + int xexp = convert_int(ax >> EXPSHIFTBITS_DP64);<br class="">
> + int xexp1 = 11 - (int) clz(ax & MANTBITS_DP64);<br class="">
> + xexp1 = xexp < 1 ? xexp1 : xexp;<br class="">
> +<br class="">
> + ulong uy = as_ulong(y);<br class="">
> + ulong ay = uy & ~SIGNBIT_DP64;<br class="">
> + double dy = as_double(ay);<br class="">
> + int yexp = convert_int(ay >> EXPSHIFTBITS_DP64);<br class="">
> + int yexp1 = 11 - (int) clz(ay & MANTBITS_DP64);<br class="">
> + yexp1 = yexp < 1 ? yexp1 : yexp;<br class="">
> +<br class="">
> + int qsgn = ((ux ^ uy) & SIGNBIT_DP64) == 0UL ? 1 : -1;<br class="">
> +<br class="">
> + // First assume |x| > |y|<br class="">
> +<br class="">
> + // Set ntimes to the number of times we need to do a<br class="">
> + // partial remainder. If the exponent of x is an exact multiple<br class="">
> + // of 53 larger than the exponent of y, and the mantissa of x is<br class="">
> + // less than the mantissa of y, ntimes will be one too large<br class="">
> + // but it doesn't matter - it just means that we'll go round<br class="">
> + // the loop below one extra time.<br class="">
> + int ntimes = max(0, (xexp1 - yexp1) / 53);<br class="">
> + double w = _clc_remainder_ldexp(dy, ntimes * 53);<br class="">
> + w = ntimes == 0 ? dy : w;<br class="">
> + double scale = ntimes == 0 ? 1.0 : 0x1.0p-53;<br class="">
> +<br class="">
> + // Each time round the loop we compute a partial remainder.<br class="">
> + // This is done by subtracting a large multiple of w<br class="">
> + // from x each time, where w is a scaled up version of y.<br class="">
> + // The subtraction must be performed exactly in quad<br class="">
> + // precision, though the result at each stage can<br class="">
> + // fit exactly in a double precision number.<br class="">
> + int i;<br class="">
> + double t, v, p, pp;<br class="">
> +<br class="">
> + for (i = 0; i < ntimes; i++) {<br class="">
> + // Compute integral multiplier<br class="">
> + t = trunc(dx / w);<br class="">
> +<br class="">
> + // Compute w * t in quad precision<br class="">
> + p = w * t;<br class="">
> + pp = fma(w, t, -p);<br class="">
> +<br class="">
> + // Subtract w * t from dx<br class="">
> + v = dx - p;<br class="">
> + dx = v + (((dx - v) - p) - pp);<br class="">
> +<br class="">
> + // If t was one too large, dx will be negative. Add back one w.<br class="">
> + dx += dx < 0.0 ? w : 0.0;<br class="">
> +<br class="">
> + // Scale w down by 2^(-53) for the next iteration<br class="">
> + w *= scale;<br class="">
> + }<br class="">
> +<br class="">
> + // One more time<br class="">
> + // Variable todd says whether the integer t is odd or not<br class="">
> + t = floor(dx / w);<br class="">
> + long lt = (long) t;<br class="">
> + int todd = lt & 1;<br class="">
> +<br class="">
> + p = w * t;<br class="">
> + pp = fma(w, t, -p);<br class="">
> + v = dx - p;<br class="">
> + dx = v + (((dx - v) - p) - pp);<br class="">
> + i = dx < 0.0;<br class="">
> + todd ^= i;<br class="">
> + dx += i ? w : 0.0;<br class="">
> +<br class="">
> + // At this point, dx lies in the range [0,dy)<br class="">
> +<br class="">
> + // For the fmod function, we're done apart from setting the correct sign.<br class="">
> + //<br class="">
> + // For the remainder function, we need to adjust dx<br class="">
> + // so that it lies in the range (-y/2, y/2] by carefully<br class="">
> + // subtracting w (== dy == y) if necessary. The rigmarole<br class="">
> + // with todd is to get the correct sign of the result<br class="">
> + // when x/y lies exactly half way between two integers,<br class="">
> + // when we need to choose the even integer.<br class="">
> +<br class="">
> + int al = (2.0 * dx > w) | (todd & (2.0 * dx == w));<br class="">
> + double dxl = dx - (al ? w : 0.0);<br class="">
> +<br class="">
> + int ag = (dx > 0.5 * w) | (todd & (dx == 0.5 * w));<br class="">
> + double dxg = dx - (ag ? w : 0.0);<br class="">
> +<br class="">
> + dx = dy < 0x1.0p+1022 ? dxl : dxg;<br class="">
> +<br class="">
> + double ret = as_double(xsgn ^ as_ulong(dx));<br class="">
> + dx = as_double(ax);<br class="">
> +<br class="">
> + // Now handle |x| == |y|<br class="">
> + int c = dx == dy;<br class="">
> + t = as_double(xsgn);<br class="">
> + ret = c ? t : ret;<br class="">
> +<br class="">
> + // Next, handle |x| < |y|<br class="">
> + c = dx < dy;<br class="">
> + ret = c ? x : ret;<br class="">
> +<br class="">
> + c &= (yexp < 1023 & 2.0 * dx > dy) | (dx > 0.5 * dy);<br class="">
> + // we could use a conversion here instead since qsgn = +-1<br class="">
> + p = qsgn == 1 ? -1.0 : 1.0;<br class="">
> + t = fma(y, p, x);<br class="">
> + ret = c ? t : ret;<br class="">
> +<br class="">
> + // We don't need anything special for |x| == 0<br class="">
> +<br class="">
> + // |y| is 0<br class="">
> + c = dy == 0.0;<br class="">
> + ret = c ? as_double(QNANBITPATT_DP64) : ret;<br class="">
> +<br class="">
> + // y is +-Inf, NaN<br class="">
> + c = yexp > BIASEDEMAX_DP64;<br class="">
> + t = y == y ? x : y;<br class="">
> + ret = c ? t : ret;<br class="">
> +<br class="">
> + // x is +=Inf, NaN<br class="">
> + c = xexp > BIASEDEMAX_DP64;<br class="">
> + ret = c ? as_double(QNANBITPATT_DP64) : ret;<br class="">
> +<br class="">
> + return ret;<br class="">
> +}<br class="">
> +<br class="">
> +_CLC_BINARY_VECTORIZE(_CLC_<wbr class="">OVERLOAD _CLC_DEF, double, remainder, double, double);<br class="">
> +#endif<br class="">
> \ No newline at end of file<br class="">
> --<br class="">
> 2.9.3<br class="">
><br class="">
> ______________________________<wbr class="">_________________<br class="">
> Libclc-dev mailing list<br class="">
> <a href="mailto:Libclc-dev@lists.llvm.org" class="">Libclc-dev@lists.llvm.org</a><br class="">
> <a href="http://lists.llvm.org/cgi-bin/mailman/listinfo/libclc-dev" rel="noreferrer" target="_blank" class="">http://lists.llvm.org/cgi-bin/<wbr class="">mailman/listinfo/libclc-dev</a><br class="">
<br class="">
<br class=""></blockquote></div><br class=""></div></div>
</div></blockquote></div><br class=""></body></html>