[PATCH 1/4] pow: Port from amd_builtins

Passes piglit on turks and carrizo
fp64 passes CTS on carrizo

Signed-off-by: Jan Vesely <jan.vesely@rutgers.edu>

Passes piglit on turks and carrizo
fp64 passes CTS on carrizo

Signed-off-by: Jan Vesely <jan.vesely@rutgers.edu>

Passes piglit on turks and carrizo
fp64 passes cts on carrizo

Signed-off-by: Jan Vesely <jan.vesely@rutgers.edu>

Passes piglit on turks and carrizo
fp64 passes ctx on carrizo

Signed-off-by: Jan Vesely <jan.vesely@rutgers.edu>

Hi Jan,

This LGTM. Two small comments below.

Passes piglit on turks and carrizo
fp64 passes CTS on carrizo

Signed-off-by: Jan Vesely <jan.vesely@rutgers.edu>
---
fp32 cts failure on carrizo is denormal related

Could you clarify this? I see some ifdefs that give the impression that they should handle this.

Jeroen

generic/include/clc/math/pown.h | 27 +--
generic/include/clc/math/pown.inc | 1 +
generic/include/math/clc_pown.h | 3 +
generic/include/math/clc_pown.inc | 1 +
generic/lib/SOURCES | 1 +
generic/lib/math/clc_pown.cl | 378 ++++++++++++++++++++++++++++++++++++++
generic/lib/math/pown.cl | 10 +-
generic/lib/math/pown.inc | 3 +
8 files changed, 393 insertions(+), 31 deletions(-)
create mode 100644 generic/include/clc/math/pown.inc
create mode 100644 generic/include/math/clc_pown.h
create mode 100644 generic/include/math/clc_pown.inc
create mode 100644 generic/lib/math/clc_pown.cl
create mode 100644 generic/lib/math/pown.inc

diff --git a/generic/include/clc/math/pown.h b/generic/include/clc/math/pown.h
index bdbf50c..2ed8a6e 100644
--- a/generic/include/clc/math/pown.h
+++ b/generic/include/clc/math/pown.h
@@ -1,24 +1,3 @@
-#define _CLC_POWN_INTRINSIC "llvm.powi"
-
-#define _CLC_POWN_DECL(GENTYPE, INTTYPE) \
- _CLC_OVERLOAD _CLC_DECL GENTYPE pown(GENTYPE x, INTTYPE y);
-
-#define _CLC_VECTOR_POWN_DECL(GENTYPE, INTTYPE) \
- _CLC_POWN_DECL(GENTYPE##2, INTTYPE##2) \
- _CLC_POWN_DECL(GENTYPE##3, INTTYPE##3) \
- _CLC_POWN_DECL(GENTYPE##4, INTTYPE##4) \
- _CLC_POWN_DECL(GENTYPE##8, INTTYPE##8) \
- _CLC_POWN_DECL(GENTYPE##16, INTTYPE##16)
-
-_CLC_OVERLOAD float pown(float x, int y) __asm(_CLC_POWN_INTRINSIC ".f32");
-
-_CLC_VECTOR_POWN_DECL(float, int)
-
-#ifdef cl_khr_fp64
-_CLC_OVERLOAD double pown(double x, int y) __asm(_CLC_POWN_INTRINSIC ".f64");
-_CLC_VECTOR_POWN_DECL(double, int)
-#endif
-
-#undef _CLC_POWN_INTRINSIC
-#undef _CLC_POWN_DECL
-#undef _CLC_VECTOR_POWN_DECL
+#define __CLC_BODY <clc/math/pown.inc>
+#include <clc/math/gentype.inc>
+#undef __CLC_BODY
diff --git a/generic/include/clc/math/pown.inc b/generic/include/clc/math/pown.inc
new file mode 100644
index 0000000..cf0be4c
--- /dev/null
+++ b/generic/include/clc/math/pown.inc
@@ -0,0 +1 @@
+_CLC_OVERLOAD _CLC_DECL __CLC_GENTYPE pown(__CLC_GENTYPE a, __CLC_INTN b);
diff --git a/generic/include/math/clc_pown.h b/generic/include/math/clc_pown.h
new file mode 100644
index 0000000..222c1dd
--- /dev/null
+++ b/generic/include/math/clc_pown.h
@@ -0,0 +1,3 @@
+#define __CLC_BODY <math/clc_pown.inc>
+#include <clc/math/gentype.inc>
+#undef __CLC_BODY
diff --git a/generic/include/math/clc_pown.inc b/generic/include/math/clc_pown.inc
new file mode 100644
index 0000000..168cd8e
--- /dev/null
+++ b/generic/include/math/clc_pown.inc
@@ -0,0 +1 @@
+_CLC_OVERLOAD _CLC_DECL __CLC_GENTYPE __clc_pown(__CLC_GENTYPE a, __CLC_INTN b);
diff --git a/generic/lib/SOURCES b/generic/lib/SOURCES
index a2d8c70..e4bd381 100644
--- a/generic/lib/SOURCES
+++ b/generic/lib/SOURCES
@@ -153,6 +153,7 @@ math/clc_nextafter.cl
math/nextafter.cl
math/clc_pow.cl
math/pow.cl
+math/clc_pown.cl
math/pown.cl
math/sin.cl
math/sincos.cl
diff --git a/generic/lib/math/clc_pown.cl b/generic/lib/math/clc_pown.cl
new file mode 100644
index 0000000..12c77a2
--- /dev/null
+++ b/generic/lib/math/clc_pown.cl
@@ -0,0 +1,378 @@
+/*
+ * 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.h"
+#include "tables.h"
+#include "../clcmacro.h"
+
+// compute pow using log and exp
+// x^y = exp(y * log(x))
+//
+// we take care not to lose precision in the intermediate steps
+//
+// When computing log, calculate it in splits,
+//
+// r = f * (p_invead + p_inv_tail)
+// r = rh + rt
+//
+// calculate log polynomial using r, in end addition, do
+// poly = poly + ((rh-r) + rt)
+//
+// lth = -r
+// ltt = ((xexp * log2_t) - poly) + logT
+// lt = lth + ltt
+//
+// lh = (xexp * log2_h) + logH
+// l = lh + lt
+//
+// Calculate final log answer as gh and gt,
+// gh = l & higher-half bits
+// gt = (((ltt - (lt - lth)) + ((lh - l) + lt)) + (l - gh))
+//
+// yh = y & higher-half bits
+// yt = y - yh
+//
+// Before entering computation of exp,
+// vs = ((yt*gt + yt*gh) + yh*gt)
+// v = vs + yh*gh
+// vt = ((yh*gh - v) + vs)
+//
+// In calculation of exp, add vt to r that is used for poly
+// At the end of exp, do
+// ((((expT * poly) + expT) + expH*poly) + expH)
+
+_CLC_DEF _CLC_OVERLOAD float __clc_pown(float x, int ny)
+{
+ float y = (float)ny;
+
+ int ix = as_int(x);
+ int ax = ix & EXSIGNBIT_SP32;
+ int xpos = ix == ax;
+
+ int iy = as_int(y);
+ int ay = iy & EXSIGNBIT_SP32;
+ int ypos = iy == ay;
+
+ // Extra precise log calculation
+ // First handle case that x is close to 1
+ float r = 1.0f - as_float(ax);
+ int near1 = fabs(r) < 0x1.0p-4f;
+ float r2 = r*r;
+
+ // Coefficients are just 1/3, 1/4, 1/5 and 1/6
+ float poly = mad(r,
+ mad(r,
+ mad(r,
+ mad(r, 0x1.24924ap-3f, 0x1.555556p-3f),
+ 0x1.99999ap-3f),
+ 0x1.000000p-2f),
+ 0x1.555556p-2f);
+
+ poly *= r2*r;
+
+ float lth_near1 = -r2 * 0.5f;
+ float ltt_near1 = -poly;
+ float lt_near1 = lth_near1 + ltt_near1;
+ float lh_near1 = -r;
+ float l_near1 = lh_near1 + lt_near1;
+
+ // Computations for x not near 1
+ int m = (int)(ax >> EXPSHIFTBITS_SP32) - EXPBIAS_SP32;
+ float mf = (float)m;
+ int ixs = as_int(as_float(ax | 0x3f800000) - 1.0f);
+ float mfs = (float)((ixs >> EXPSHIFTBITS_SP32) - 253);
+ int c = m == -127;
+ int ixn = c ? ixs : ax;
+ float mfn = c ? mfs : mf;
+
+ int indx = (ixn & 0x007f0000) + ((ixn & 0x00008000) << 1);
+
+ // F - Y
+ float f = as_float(0x3f000000 | indx) - as_float(0x3f000000 | (ixn & MANTBITS_SP32));
+
+ indx = indx >> 16;
+ float2 tv = USE_TABLE(log_inv_tbl_ep, indx);
+ float rh = f * tv.s0;
+ float rt = f * tv.s1;
+ r = rh + rt;
+
+ poly = mad(r, mad(r, 0x1.0p-2f, 0x1.555556p-2f), 0x1.0p-1f) * (r*r);
+ poly += (rh - r) + rt;
+
+ const float LOG2_HEAD = 0x1.62e000p-1f; // 0.693115234
+ const float LOG2_TAIL = 0x1.0bfbe8p-15f; // 0.0000319461833
+ tv = USE_TABLE(loge_tbl, indx);
+ float lth = -r;
+ float ltt = mad(mfn, LOG2_TAIL, -poly) + tv.s1;
+ float lt = lth + ltt;
+ float lh = mad(mfn, LOG2_HEAD, tv.s0);
+ float l = lh + lt;
+
+ // Select near 1 or not
+ lth = near1 ? lth_near1 : lth;
+ ltt = near1 ? ltt_near1 : ltt;
+ lt = near1 ? lt_near1 : lt;
+ lh = near1 ? lh_near1 : lh;
+ l = near1 ? l_near1 : l;
+
+ float gh = as_float(as_int(l) & 0xfffff000);
+ float gt = ((ltt - (lt - lth)) + ((lh - l) + lt)) + (l - gh);
+
+ float yh = as_float(iy & 0xfffff000);
+
+ float yt = (float)(ny - (int)yh);
+
+ float ylogx_s = mad(gt, yh, mad(gh, yt, yt*gt));
+ float ylogx = mad(yh, gh, ylogx_s);
+ float ylogx_t = mad(yh, gh, -ylogx) + ylogx_s;
+
+ // Extra precise exp of ylogx
+ const float R_64_BY_LOG2 = 0x1.715476p+6f; // 64/log2 : 92.332482616893657
+ int n = convert_int(ylogx * R_64_BY_LOG2);
+ float nf = (float) n;
+
+ int j = n & 0x3f;
+ m = n >> 6;
+ int m2 = m << EXPSHIFTBITS_SP32;
+
+ const float R_LOG2_BY_64_LD = 0x1.620000p-7f; // log2/64 lead: 0.0108032227
+ const float R_LOG2_BY_64_TL = 0x1.c85fdep-16f; // log2/64 tail: 0.0000272020388
+ r = mad(nf, -R_LOG2_BY_64_TL, mad(nf, -R_LOG2_BY_64_LD, ylogx)) + ylogx_t;
+
+ // Truncated Taylor series for e^r
+ poly = mad(mad(mad(r, 0x1.555556p-5f, 0x1.555556p-3f), r, 0x1.000000p-1f), r*r, r);
+
+ tv = USE_TABLE(exp_tbl_ep, j);
+
+ float expylogx = mad(tv.s0, poly, mad(tv.s1, poly, tv.s1)) + tv.s0;
+ #if !defined(SUBNORMALS_SUPPORTED)
+ int explg = ((as_uint(expylogx) & EXPBITS_SP32 >> 23) - 127);
+ m = (23-(m + 149)) == 0 ? 1: m;
+ uint mantissa = ((as_uint(expylogx) & MANTBITS_SP32)|IMPBIT_SP32) >> (23-(m + 149));
+ float sexpylogx = as_float(mantissa);
+ #else
+ float sexpylogx = expylogx * as_float(0x1 << (m + 149));
+ #endif
+
+
+ float texpylogx = as_float(as_int(expylogx) + m2);
+ expylogx = m < -125 ? sexpylogx : texpylogx;
+
+ // Result is +-Inf if (ylogx + ylogx_t) > 128*log2
+ expylogx = ((ylogx > 0x1.62e430p+6f) | (ylogx == 0x1.62e430p+6f & ylogx_t > -0x1.05c610p-22f)) ? as_float(PINFBITPATT_SP32) : expylogx;
+
+ // Result is 0 if ylogx < -149*log2
+ expylogx = ylogx < -0x1.9d1da0p+6f ? 0.0f : expylogx;
+
+ // Classify y:
+ // inty = 0 means not an integer.
+ // inty = 1 means odd integer.
+ // inty = 2 means even integer.
+
+ int inty = 2 - (ny & 1);
+
+ float signval = as_float((as_uint(expylogx) ^ SIGNBIT_SP32));
+ expylogx = ((inty == 1) & !xpos) ? signval : expylogx;
+ int ret = as_int(expylogx);
+
+ // Corner case handling
+ int xinf = xpos ? PINFBITPATT_SP32 : NINFBITPATT_SP32;
+ ret = ((ax == 0) & !ypos & (inty == 1)) ? xinf : ret;
+ ret = ((ax == 0) & !ypos & (inty == 2)) ? PINFBITPATT_SP32 : ret;
+ ret = ((ax == 0) & ypos & (inty == 2)) ? 0 : ret;
+ int xzero = !xpos ? 0x80000000 : 0L;
+ ret = ((ax == 0) & ypos & (inty == 1)) ? xzero : ret;
+ ret = ((ix == NINFBITPATT_SP32) & !ypos & (inty == 1)) ? 0x80000000 : ret;
+ ret = ((ix == NINFBITPATT_SP32) & !ypos & (inty != 1)) ? 0 : ret;
+ ret = ((ix == NINFBITPATT_SP32) & ypos & (inty == 1)) ? NINFBITPATT_SP32 : ret;
+ ret = ((ix == NINFBITPATT_SP32) & ypos & (inty != 1)) ? PINFBITPATT_SP32 : ret;
+ ret = ((ix == PINFBITPATT_SP32) & !ypos) ? 0 : ret;
+ ret = ((ix == PINFBITPATT_SP32) & ypos) ? PINFBITPATT_SP32 : ret;
+ ret = ax > PINFBITPATT_SP32 ? ix : ret;
+ ret = ny == 0 ? 0x3f800000 : ret;
+
+ return as_float(ret);
+}
+_CLC_BINARY_VECTORIZE(_CLC_DEF _CLC_OVERLOAD, float, __clc_pown, float, int)
+
+#ifdef cl_khr_fp64
+_CLC_DEF _CLC_OVERLOAD double __clc_pown(double x, int ny)
+{
+ const double real_log2_tail = 5.76999904754328540596e-08;
+ const double real_log2_lead = 6.93147122859954833984e-01;
+
+ double y = (double) ny;
+
+ long ux = as_long(x);
+ long ax = ux & (~SIGNBIT_DP64);
+ int xpos = ax == ux;
+
+ long uy = as_long(y);
+ long ay = uy & (~SIGNBIT_DP64);
+ int ypos = ay == uy;
+
+ // Extended precision log
+ double v, vt;
+ {
+ int exp = (int)(ax >> 52) - 1023;
+ int mask_exp_1023 = exp == -1023;
+ double xexp = (double) exp;
+ long mantissa = ax & 0x000FFFFFFFFFFFFFL;
+
+ long temp_ux = as_long(as_double(0x3ff0000000000000L | mantissa) - 1.0);
+ exp = ((temp_ux & 0x7FF0000000000000L) >> 52) - 2045;
+ double xexp1 = (double) exp;
+ long mantissa1 = temp_ux & 0x000FFFFFFFFFFFFFL;
+
+ xexp = mask_exp_1023 ? xexp1 : xexp;
+ mantissa = mask_exp_1023 ? mantissa1 : mantissa;
+
+ long rax = (mantissa & 0x000ff00000000000) + ((mantissa & 0x0000080000000000) << 1);
+ int index = rax >> 44;
+
+ double F = as_double(rax | 0x3FE0000000000000L);
+ double Y = as_double(mantissa | 0x3FE0000000000000L);
+ double f = F - Y;
+ double2 tv = USE_TABLE(log_f_inv_tbl, index);
+ double log_h = tv.s0;
+ double log_t = tv.s1;
+ double f_inv = (log_h + log_t) * f;
+ double r1 = as_double(as_long(f_inv) & 0xfffffffff8000000L);
+ double r2 = fma(-F, r1, f) * (log_h + log_t);
+ double r = r1 + r2;
+
+ double poly = fma(r,
+ fma(r,
+ fma(r,
+ fma(r, 1.0/7.0, 1.0/6.0),
+ 1.0/5.0),
+ 1.0/4.0),
+ 1.0/3.0);
+ poly = poly * r * r * r;
+
+ double hr1r1 = 0.5*r1*r1;
+ double poly0h = r1 + hr1r1;
+ double poly0t = r1 - poly0h + hr1r1;
+ poly = fma(r1, r2, fma(0.5*r2, r2, poly)) + r2 + poly0t;

There seems something wrong with the indentation here. The same holds for the other two patches.

Jeroen

Was there a patch 1 for this series? Patch 2 doesn't apply cleanly as-
is on top of current HEAD/trunk/tip/whatever.

--Aaron

Was there a patch 1 for this series? Patch 2 doesn't apply cleanly as-
is on top of current HEAD/trunk/tip/whatever.

there was, the ML complained about it's size (it adds few tables to
tables.cl), but I thought it went through.
the whole series can be find here:
https://github.com/jvesely/libclc/tree/pow

Jan

Hi Jan,

This LGTM. Two small comments below.

>
> Passes piglit on turks and carrizo
> fp64 passes CTS on carrizo
>
> Signed-off-by: Jan Vesely <jan.vesely@rutgers.edu>
> ---
> fp32 cts failure on carrizo is denormal related

Could you clarify this? I see some ifdefs that give the impression that they should handle this.

thanks for bringing this up. Those should be converted to runtime
checks. I'll send v2.

I see errors like this:
ERROR: rootn: -2047.000000 ulp error at {-0x1.ffcp-139, 1}: *-0x1.ffcp-
139 vs. -0x1.ffcp-138 (0x80000ffe) at index: 137

ERROR: powr: 8106.990234 ulp error at {0x1.e70f2cp-18 (0x36f38796),
0x1.fde54cp+2 (0x40fef2a6)}: *0x1.fabp-137 vs. 0x1.fabp-136
(0x00003f56) at index: 64996

ERROR: pown: -2047.000000 ulp error at {-0x1.ffcp-139, 1}: *-0x1.ffcp-
139 vs. -0x1.ffcp-138 (0x80000ffe) at index: 137

ERROR: pow: 8106.990234 ulp error at {0x1.e70f2cp-18 (0x36f38796),
0x1.fde54cp+2 (0x40fef2a6)}: *0x1.fabp-137 vs. 0x1.fabp-136
(0x00003f56) at index: 64996

rootn and pown should probably be specialcased (y = 1), although the
specs do not require it, but there might be more similar errors.
I'm not sure what to do about pow/powr it looks like off by one error
in denormal path, but I think that would be better addressed in a
separate patch.

Jan

PS: I'll fix the formatting in v2 as well.

Hi Jan,

This LGTM. Two small comments below.

Passes piglit on turks and carrizo
fp64 passes CTS on carrizo

Signed-off-by: Jan Vesely <jan.vesely@rutgers.edu>

fp32 cts failure on carrizo is denormal related

Could you clarify this? I see some ifdefs that give the impression that they should handle this.

thanks for bringing this up. Those should be converted to runtime
checks. I’ll send v2.

I see errors like this:
ERROR: rootn: -2047.000000 ulp error at {-0x1.ffcp-139, 1}: *-0x1.ffcp-
139 vs. -0x1.ffcp-138 (0x80000ffe) at index: 137

ERROR: powr: 8106.990234 ulp error at {0x1.e70f2cp-18 (0x36f38796),
0x1.fde54cp+2 (0x40fef2a6)}: *0x1.fabp-137 vs. 0x1.fabp-136
(0x00003f56) at index: 64996

ERROR: pown: -2047.000000 ulp error at {-0x1.ffcp-139, 1}: *-0x1.ffcp-
139 vs. -0x1.ffcp-138 (0x80000ffe) at index: 137

ERROR: pow: 8106.990234 ulp error at {0x1.e70f2cp-18 (0x36f38796),
0x1.fde54cp+2 (0x40fef2a6)}: *0x1.fabp-137 vs. 0x1.fabp-136
(0x00003f56) at index: 64996

rootn and pown should probably be specialcased (y = 1), although the
specs do not require it, but there might be more similar errors.
I’m not sure what to do about pow/powr it looks like off by one error
in denormal path, but I think that would be better addressed in a
separate patch.

Agreed that this is better addressed separately. I doubt if special-casing is
the solution here. It somehow seems like they all suffer from a similar
kind off-by-one error in the exponent.

Jeroen