libclc/generic/lib/math/tgamma.cl

   1 /*
   2  * Copyright (c) 2016 Aaron Watry
   3  * Copyright (c) 2014 Advanced Micro Devices, Inc.
   4  *
   5  * Permission is hereby granted, free of charge, to any person obtaining a copy
   6  * of this software and associated documentation files (the "Software"), to deal
   7  * in the Software without restriction, including without limitation the rights
   8  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
   9  * copies of the Software, and to permit persons to whom the Software is
  10  * furnished to do so, subject to the following conditions:
  11  *
  12  * The above copyright notice and this permission notice shall be included in
  13  * all copies or substantial portions of the Software.
  14  *
  15  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  16  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  17  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  18  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  19  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  20  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  21  * THE SOFTWARE.
  22  */
  23
  24 #include <clc/clc.h>
  25 #include <clc/clcmacro.h>
  26
  27 #include "math.h"
  28
  29 _CLC_OVERLOAD _CLC_DEF float tgamma(float x) {
  30     const float pi = 3.1415926535897932384626433832795f;
  31     float ax = fabs(x);
  32     float lg = lgamma(ax);
  33     float g = exp(lg);
  34
  35     if (x < 0.0f) {
  36         float z = sinpi(x);
  37         g = g * ax * z;
  38         g = pi / g;
  39         g = g == 0 ? as_float(PINFBITPATT_SP32) : g;
  40         g = z == 0 ? as_float(QNANBITPATT_SP32) : g;
  41     }
  42
  43     return g;
  44 }
  45
  46 _CLC_UNARY_VECTORIZE(_CLC_OVERLOAD _CLC_DEF, float, tgamma, float);
  47
  48 #ifdef cl_khr_fp64
  49
  50 #pragma OPENCL EXTENSION cl_khr_fp64 : enable
  51
  52 _CLC_OVERLOAD _CLC_DEF double tgamma(double x) {
  53     const double pi = 3.1415926535897932384626433832795;
  54     double ax = fabs(x);
  55     double lg = lgamma(ax);
  56     double g = exp(lg);
  57
  58     if (x < 0.0) {
  59         double z = sinpi(x);
  60         g = g * ax * z;
  61         g = pi / g;
  62         g = g == 0 ? as_double(PINFBITPATT_DP64) : g;
  63         g = z == 0 ? as_double(QNANBITPATT_DP64) : g;
  64     }
  65
  66     return g;
  67 }
  68
  69 _CLC_UNARY_VECTORIZE(_CLC_OVERLOAD _CLC_DEF, double, tgamma, double);
  70
  71 #endif