[Libclc-dev] [PATCH 2/4] pown: Port from amd_builtins

Jeroen Ketema via Libclc-dev libclc-dev at lists.llvm.org
Wed Jan 17 12:59:23 PST 2018



> On 17 Jan 2018, at 17:27, Jan Vesely <jan.vesely at rutgers.edu> wrote:
> 
> On Tue, 2018-01-16 at 20:32 +0100, Jeroen Ketema via Libclc-dev wrote:
>> Hi Jan,
>> 
>> This LGTM. Two small comments below.
>> 
>>> On 16 Jan 2018, at 01:47, Jan Vesely via Libclc-dev <libclc-dev at lists.llvm.org> wrote:
>>> 
>>> Passes piglit on turks and carrizo
>>> fp64 passes CTS on carrizo
>>> 
>>> Signed-off-by: Jan Vesely <jan.vesely at 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

> 
> Jan
> 
> PS: I'll fix the formatting in v2 as well.
> 
>> 
>> 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
>> 
>>> +
>>> +        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;
>>> +
>>> +        int mask_2_24 = ay > 0x4170000000000000; // 2^24
>>> +        int nyh = convert_int(y_head);
>>> +        int nyt = ny - nyh;
>>> +        double y_tail1 = (double)nyt;
>>> +        y_tail = mask_2_24 ? y_tail1 : y_tail;
>>> +
>>> +        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
>>> +    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 ? 0x8000000000000000L : 0L;
>>> +    ret = ((ax == 0L) & ypos & (inty == 1)) ? xzero : ret;
>>> +    ret = ((ux == NINFBITPATT_DP64) & !ypos & (inty == 1)) ? 0x8000000000000000L : ret;
>>> +    ret = ((ux == NINFBITPATT_DP64) & !ypos & (inty != 1)) ? 0L : ret;
>>> +    ret = ((ux == NINFBITPATT_DP64) & ypos & (inty == 1)) ? NINFBITPATT_DP64 : ret;
>>> +    ret = ((ux == NINFBITPATT_DP64) & ypos & (inty != 1)) ? PINFBITPATT_DP64 : 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 ? 0x3ff0000000000000L : ret;
>>> +
>>> +    return as_double(ret);
>>> +}
>>> +_CLC_BINARY_VECTORIZE(_CLC_DEF _CLC_OVERLOAD, double, __clc_pown, double, int)
>>> +#endif
>>> diff --git a/generic/lib/math/pown.cl b/generic/lib/math/pown.cl
>>> index f3b27d4..7652cd5 100644
>>> --- a/generic/lib/math/pown.cl
>>> +++ b/generic/lib/math/pown.cl
>>> @@ -1,10 +1,6 @@
>>> #include <clc/clc.h>
>>> -#include "../clcmacro.h"
>>> 
>>> -_CLC_BINARY_VECTORIZE(_CLC_OVERLOAD _CLC_DEF, float, pown, float, int)
>>> +#include <math/clc_pown.h>
>>> 
>>> -#ifdef cl_khr_fp64
>>> -#pragma OPENCL EXTENSION cl_khr_fp64 : enable
>>> -
>>> -_CLC_BINARY_VECTORIZE(_CLC_OVERLOAD _CLC_DEF, double, pown, double, int)
>>> -#endif
>>> +#define __CLC_BODY <pown.inc>
>>> +#include <clc/math/gentype.inc>
>>> diff --git a/generic/lib/math/pown.inc b/generic/lib/math/pown.inc
>>> new file mode 100644
>>> index 0000000..84729d9
>>> --- /dev/null
>>> +++ b/generic/lib/math/pown.inc
>>> @@ -0,0 +1,3 @@
>>> +_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE pown(__CLC_GENTYPE x, __CLC_INTN y) {
>>> +  return __clc_pown(x, y);
>>> +}
>>> -- 
>>> 2.14.3
>>> 
>>> _______________________________________________
>>> Libclc-dev mailing list
>>> Libclc-dev at lists.llvm.org
>>> http://lists.llvm.org/cgi-bin/mailman/listinfo/libclc-dev
>> 
>> _______________________________________________
>> Libclc-dev mailing list
>> Libclc-dev at lists.llvm.org <mailto:Libclc-dev at lists.llvm.org>
>> http://lists.llvm.org/cgi-bin/mailman/listinfo/libclc-dev <http://lists.llvm.org/cgi-bin/mailman/listinfo/libclc-dev>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.llvm.org/pipermail/libclc-dev/attachments/20180117/c5fae357/attachment-0001.html>


More information about the Libclc-dev mailing list