Author: jvesely Date: Wed Jan 17 13:22:14 2018 New Revision: 322763 URL: http://llvm.org/viewvc/llvm-project?rev=322763&view=rev Log: rootn: Port from amd_builtins
Passes piglit on turks and carrizo fp64 passes ctx on carrizo v2: fix formatting check fp32 denormal support at runtime Reviewer: Jeroen Ketema <j.ket...@xs4all.nl> Signed-off-by: Jan Vesely <jan.ves...@rutgers.edu> Added: libclc/trunk/generic/include/clc/math/rootn.h libclc/trunk/generic/include/clc/math/rootn.inc libclc/trunk/generic/include/math/clc_rootn.h libclc/trunk/generic/include/math/clc_rootn.inc libclc/trunk/generic/lib/math/clc_rootn.cl libclc/trunk/generic/lib/math/rootn.cl libclc/trunk/generic/lib/math/rootn.inc Modified: libclc/trunk/generic/include/clc/clc.h libclc/trunk/generic/lib/SOURCES Modified: libclc/trunk/generic/include/clc/clc.h URL: http://llvm.org/viewvc/llvm-project/libclc/trunk/generic/include/clc/clc.h?rev=322763&r1=322762&r2=322763&view=diff ============================================================================== --- libclc/trunk/generic/include/clc/clc.h (original) +++ libclc/trunk/generic/include/clc/clc.h Wed Jan 17 13:22:14 2018 @@ -91,6 +91,7 @@ #include <clc/math/pown.h> #include <clc/math/powr.h> #include <clc/math/rint.h> +#include <clc/math/rootn.h> #include <clc/math/round.h> #include <clc/math/sin.h> #include <clc/math/sincos.h> Added: libclc/trunk/generic/include/clc/math/rootn.h URL: http://llvm.org/viewvc/llvm-project/libclc/trunk/generic/include/clc/math/rootn.h?rev=322763&view=auto ============================================================================== --- libclc/trunk/generic/include/clc/math/rootn.h (added) +++ libclc/trunk/generic/include/clc/math/rootn.h Wed Jan 17 13:22:14 2018 @@ -0,0 +1,3 @@ +#define __CLC_BODY <clc/math/rootn.inc> +#include <clc/math/gentype.inc> +#undef __CLC_BODY Added: libclc/trunk/generic/include/clc/math/rootn.inc URL: http://llvm.org/viewvc/llvm-project/libclc/trunk/generic/include/clc/math/rootn.inc?rev=322763&view=auto ============================================================================== --- libclc/trunk/generic/include/clc/math/rootn.inc (added) +++ libclc/trunk/generic/include/clc/math/rootn.inc Wed Jan 17 13:22:14 2018 @@ -0,0 +1 @@ +_CLC_OVERLOAD _CLC_DECL __CLC_GENTYPE rootn(__CLC_GENTYPE a, __CLC_INTN b); Added: libclc/trunk/generic/include/math/clc_rootn.h URL: http://llvm.org/viewvc/llvm-project/libclc/trunk/generic/include/math/clc_rootn.h?rev=322763&view=auto ============================================================================== --- libclc/trunk/generic/include/math/clc_rootn.h (added) +++ libclc/trunk/generic/include/math/clc_rootn.h Wed Jan 17 13:22:14 2018 @@ -0,0 +1,3 @@ +#define __CLC_BODY <math/clc_rootn.inc> +#include <clc/math/gentype.inc> +#undef __CLC_BODY Added: libclc/trunk/generic/include/math/clc_rootn.inc URL: http://llvm.org/viewvc/llvm-project/libclc/trunk/generic/include/math/clc_rootn.inc?rev=322763&view=auto ============================================================================== --- libclc/trunk/generic/include/math/clc_rootn.inc (added) +++ libclc/trunk/generic/include/math/clc_rootn.inc Wed Jan 17 13:22:14 2018 @@ -0,0 +1 @@ +_CLC_OVERLOAD _CLC_DECL __CLC_GENTYPE __clc_rootn(__CLC_GENTYPE a, __CLC_INTN b); Modified: libclc/trunk/generic/lib/SOURCES URL: http://llvm.org/viewvc/llvm-project/libclc/trunk/generic/lib/SOURCES?rev=322763&r1=322762&r2=322763&view=diff ============================================================================== --- libclc/trunk/generic/lib/SOURCES (original) +++ libclc/trunk/generic/lib/SOURCES Wed Jan 17 13:22:14 2018 @@ -144,6 +144,8 @@ math/clc_pown.cl math/pown.cl math/clc_powr.cl math/powr.cl +math/clc_rootn.cl +math/rootn.cl math/sin.cl math/sincos.cl math/sincos_helpers.cl Added: libclc/trunk/generic/lib/math/clc_rootn.cl URL: http://llvm.org/viewvc/llvm-project/libclc/trunk/generic/lib/math/clc_rootn.cl?rev=322763&view=auto ============================================================================== --- libclc/trunk/generic/lib/math/clc_rootn.cl (added) +++ libclc/trunk/generic/lib/math/clc_rootn.cl Wed Jan 17 13:22:14 2018 @@ -0,0 +1,379 @@ +/* + * 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 "config.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_rootn(float x, int ny) +{ + float y = MATH_RECIP((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 fny = (float)ny; + float fnyh = as_float(as_int(fny) & 0xfffff000); + float fnyt = (float)(ny - (int)fnyh); + float yt = MATH_DIVIDE(mad(-fnyt, yh, mad(-fnyh, yh, 1.0f)), fny); + + 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; + float sexpylogx; + if (!__clc_fp32_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)); + sexpylogx = as_float(mantissa); + } else { + sexpylogx = expylogx * as_float(0x1 << (m + 149)); + } + + + 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 + ret = (!xpos & (inty == 2)) ? QNANBITPATT_SP32 : ret; + 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 ? 0 : 0x80000000; + ret = ((ax == 0) & ypos & (inty == 1)) ? xzero : ret; + ret = ((ix == NINFBITPATT_SP32) & ypos & (inty == 1)) ? NINFBITPATT_SP32 : ret; + ret = ((ix == NINFBITPATT_SP32) & !ypos & (inty == 1)) ? 0x80000000 : 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 ? QNANBITPATT_SP32 : ret; + + return as_float(ret); +} +_CLC_BINARY_VECTORIZE(_CLC_DEF _CLC_OVERLOAD, float, __clc_rootn, float, int) + +#ifdef cl_khr_fp64 +_CLC_DEF _CLC_OVERLOAD double __clc_rootn(double x, int ny) +{ + const double real_log2_tail = 5.76999904754328540596e-08; + const double real_log2_lead = 6.93147122859954833984e-01; + + double dny = (double)ny; + double y = 1.0 / dny; + + 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; + + tv = USE_TABLE(powlog_tbl, index); + log_h = tv.s0; + log_t = tv.s1; + + double resT_t = fma(xexp, real_log2_tail, + log_t) - poly; + double resT = resT_t - poly0h; + double resH = fma(xexp, real_log2_lead, log_h); + double resT_h = poly0h; + + double H = resT + resH; + double H_h = as_double(as_long(H) & 0xfffffffff8000000L); + double T = (resH - H + resT) + (resT_t - (resT + resT_h)) + (H - H_h); + H = H_h; + + double y_head = as_double(uy & 0xfffffffff8000000L); + double y_tail = y - y_head; + + double fnyh = as_double(as_long(dny) & 0xfffffffffff00000); + double fnyt = (double)(ny - (int)fnyh); + y_tail = fma(-fnyt, y_head, fma(-fnyh, y_head, 1.0))/ dny; + + double temp = fma(y_tail, H, fma(y_head, T, y_tail*T)); + v = fma(y_head, H, temp); + vt = fma(y_head, H, -v) + temp; + } + + // Now calculate exp of (v,vt) + + double expv; + { + const double max_exp_arg = 709.782712893384; + const double min_exp_arg = -745.1332191019411; + const double sixtyfour_by_lnof2 = 92.33248261689366; + const double lnof2_by_64_head = 0.010830424260348081; + const double lnof2_by_64_tail = -4.359010638708991e-10; + + double temp = v * sixtyfour_by_lnof2; + int n = (int)temp; + double dn = (double)n; + int j = n & 0x0000003f; + int m = n >> 6; + + double2 tv = USE_TABLE(two_to_jby64_ep_tbl, j); + double f1 = tv.s0; + double f2 = tv.s1; + double f = f1 + f2; + + double r1 = fma(dn, -lnof2_by_64_head, v); + double r2 = dn * lnof2_by_64_tail; + double r = (r1 + r2) + vt; + + double q = fma(r, + fma(r, + fma(r, + fma(r, 1.38889490863777199667e-03, 8.33336798434219616221e-03), + 4.16666666662260795726e-02), + 1.66666666665260878863e-01), + 5.00000000000000008883e-01); + q = fma(r*r, q, r); + + expv = fma(f, q, f2) + f1; + expv = ldexp(expv, m); + + expv = v > max_exp_arg ? as_double(0x7FF0000000000000L) : expv; + expv = v < min_exp_arg ? 0.0 : expv; + } + + // See whether y is an integer. + // inty = 0 means not an integer. + // inty = 1 means odd integer. + // inty = 2 means even integer. + + int inty = 2 - (ny & 1); + + expv *= ((inty == 1) & !xpos) ? -1.0 : 1.0; + + long ret = as_long(expv); + + // Now all the edge cases + ret = (!xpos & (inty == 2)) ? QNANBITPATT_DP64 : ret; + long xinf = xpos ? PINFBITPATT_DP64 : NINFBITPATT_DP64; + ret = ((ax == 0L) & !ypos & (inty == 1)) ? xinf : ret; + ret = ((ax == 0L) & !ypos & (inty == 2)) ? PINFBITPATT_DP64 : ret; + ret = ((ax == 0L) & ypos & (inty == 2)) ? 0L : ret; + long xzero = xpos ? 0L : 0x8000000000000000L; + ret = ((ax == 0L) & ypos & (inty == 1)) ? xzero : ret; + ret = ((ux == NINFBITPATT_DP64) & ypos & (inty == 1)) ? NINFBITPATT_DP64 : ret; + ret = ((ux == NINFBITPATT_DP64) & !ypos & (inty == 1)) ? 0x8000000000000000L : ret; + ret = ((ux == PINFBITPATT_DP64) & !ypos) ? 0L : ret; + ret = ((ux == PINFBITPATT_DP64) & ypos) ? PINFBITPATT_DP64 : ret; + ret = ax > PINFBITPATT_DP64 ? ux : ret; + ret = ny == 0 ? QNANBITPATT_DP64 : ret; + return as_double(ret); +} +_CLC_BINARY_VECTORIZE(_CLC_DEF _CLC_OVERLOAD, double, __clc_rootn, double, int) +#endif Added: libclc/trunk/generic/lib/math/rootn.cl URL: http://llvm.org/viewvc/llvm-project/libclc/trunk/generic/lib/math/rootn.cl?rev=322763&view=auto ============================================================================== --- libclc/trunk/generic/lib/math/rootn.cl (added) +++ libclc/trunk/generic/lib/math/rootn.cl Wed Jan 17 13:22:14 2018 @@ -0,0 +1,6 @@ +#include <clc/clc.h> + +#include <math/clc_rootn.h> + +#define __CLC_BODY <rootn.inc> +#include <clc/math/gentype.inc> Added: libclc/trunk/generic/lib/math/rootn.inc URL: http://llvm.org/viewvc/llvm-project/libclc/trunk/generic/lib/math/rootn.inc?rev=322763&view=auto ============================================================================== --- libclc/trunk/generic/lib/math/rootn.inc (added) +++ libclc/trunk/generic/lib/math/rootn.inc Wed Jan 17 13:22:14 2018 @@ -0,0 +1,3 @@ +_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE rootn(__CLC_GENTYPE x, __CLC_INTN y) { + return __clc_rootn(x, y); +} _______________________________________________ cfe-commits mailing list cfe-commits@lists.llvm.org http://lists.llvm.org/cgi-bin/mailman/listinfo/cfe-commits