libclc/generic/lib/math/cospi.cl

   1 /*
   2  * Copyright (c) 2014 Advanced Micro Devices, Inc.
   3  *
   4  * Permission is hereby granted, free of charge, to any person obtaining a copy
   5  * of this software and associated documentation files (the "Software"), to deal
   6  * in the Software without restriction, including without limitation the rights
   7  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
   8  * copies of the Software, and to permit persons to whom the Software is
   9  * furnished to do so, subject to the following conditions:
  10  *
  11  * The above copyright notice and this permission notice shall be included in
  12  * all copies or substantial portions of the Software.
  13  *
  14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  15  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  16  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  17  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  18  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  19  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  20  * THE SOFTWARE.
  21  */
  22
  23 #include <clc/clc.h>
  24
  25 #include "math.h"
  26 #include "sincos_helpers.h"
  27 #include "sincospiF_piby4.h"
  28 #include "../clcmacro.h"
  29 #ifdef cl_khr_fp64
  30 #include "sincosD_piby4.h"
  31 #endif
  32
  33 _CLC_OVERLOAD _CLC_DEF float cospi(float x)
  34 {
  35     int ix = as_int(x) & 0x7fffffff;
  36     float ax = as_float(ix);
  37     int iax = (int)ax;
  38     float r = ax - iax;
  39     int xodd = iax & 0x1 ? 0x80000000 : 0;
  40
  41     // Initialize with return for +-Inf and NaN
  42     int ir = 0x7fc00000;
  43
  44     // 2^24 <= |x| < Inf, the result is always even integer
  45     ir = ix < 0x7f800000 ? 0x3f800000 : ir;
  46
  47     // 2^23 <= |x| < 2^24, the result is always integer
  48     ir = ix < 0x4b800000 ? xodd | 0x3f800000 : ir;
  49
  50     // 0x1.0p-7 <= |x| < 2^23, result depends on which 0.25 interval
  51
  52     // r < 1.0
  53     float a = 1.0f - r;
  54     int e = 1;
  55     int s = xodd ^ 0x80000000;
  56
  57     // r <= 0.75
  58     int c = r <= 0.75f;
  59     a = c ? r - 0.5f : a;
  60     e = c ? 0 : e;
  61
  62     // r < 0.5
  63     c = r < 0.5f;
  64     a = c ? 0.5f - r : a;
  65     s = c ? xodd : s;
  66
  67     // r <= 0.25
  68     c = r <= 0.25f;
  69     a = c ? r : a;
  70     e = c ? 1 : e;
  71
  72     float2 t = __libclc__sincosf_piby4(a * M_PI_F);
  73     int jr = s ^ as_int(e ? t.hi : t.lo);
  74
  75     ir = ix < 0x4b000000 ? jr : ir;
  76
  77     return as_float(ir);
  78 }
  79
  80
  81 _CLC_UNARY_VECTORIZE(_CLC_OVERLOAD _CLC_DEF, float, cospi, float);
  82
  83 #ifdef cl_khr_fp64
  84
  85 #pragma OPENCL EXTENSION cl_khr_fp64 : enable
  86
  87 _CLC_OVERLOAD _CLC_DEF double cospi(double x) {
  88
  89     long ix = as_long(x) & 0x7fffffffffffffffL;
  90     double ax = as_double(ix);
  91     long iax = (long)ax;
  92     double r = ax - (double)iax;
  93     long xodd = iax & 0x1L ? 0x8000000000000000L : 0L;
  94
  95     // Initialize with return for +-Inf and NaN
  96     long ir = 0x7ff8000000000000L;
  97
  98     // 2^53 <= |x| < Inf, the result is always even integer
  99     ir = ix < 0x7ff0000000000000 ? 0x3ff0000000000000L : ir;
 100
 101     // 2^52 <= |x| < 2^53, the result is always integer
 102     ir = ax < 0x1.0p+53 ? xodd | 0x3ff0000000000000L : ir;
 103
 104     // 0x1.0p-7 <= |x| < 2^52, result depends on which 0.25 interval
 105
 106     // r < 1.0
 107     double a = 1.0 - r;
 108     int e = 1;
 109     long s = xodd ^ 0x8000000000000000L;
 110
 111     // r <= 0.75
 112     int c = r <= 0.75;
 113     double t = r - 0.5;
 114     a = c ? t : a;
 115     e = c ? 0 : e;
 116
 117     // r < 0.5
 118     c = r < 0.5;
 119     t = 0.5 - r;
 120     a = c ? t : a;
 121     s = c ? xodd : s;
 122
 123     // r <= 0.25
 124     c = r <= 0.25;
 125     a = c ? r : a;
 126     e = c ? 1 : e;
 127
 128     double2 sc = __libclc__sincos_piby4(a * M_PI, 0.0);
 129     long jr = s ^ as_long(e ? sc.hi : sc.lo);
 130
 131     ir = ax < 0x1.0p+52 ? jr : ir;
 132
 133     return as_double(ir);
 134 }
 135 _CLC_UNARY_VECTORIZE(_CLC_OVERLOAD _CLC_DEF, double, cospi, double);
 136 #endif