libclc/generic/lib/math/atan2.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 "tables.h"
  27 #include "../clcmacro.h"
  28
  29 _CLC_OVERLOAD _CLC_DEF float atan2(float y, float x)
  30 {
  31     const float pi = 0x1.921fb6p+1f;
  32     const float piby2 = 0x1.921fb6p+0f;
  33     const float piby4 = 0x1.921fb6p-1f;
  34     const float threepiby4 = 0x1.2d97c8p+1f;
  35
  36     float ax = fabs(x);
  37     float ay = fabs(y);
  38     float v = min(ax, ay);
  39     float u = max(ax, ay);
  40
  41     // Scale since u could be large, as in "regular" divide
  42     float s = u > 0x1.0p+96f ? 0x1.0p-32f : 1.0f;
  43     float vbyu = s * MATH_DIVIDE(v, s*u);
  44
  45     float vbyu2 = vbyu * vbyu;
  46
  47 #define USE_2_2_APPROXIMATION
  48 #if defined USE_2_2_APPROXIMATION
  49     float p = mad(vbyu2, mad(vbyu2, -0x1.7e1f78p-9f, -0x1.7d1b98p-3f), -0x1.5554d0p-2f) * vbyu2 * vbyu;
  50     float q = mad(vbyu2, mad(vbyu2, 0x1.1a714cp-2f, 0x1.287c56p+0f), 1.0f);
  51 #else
  52     float p = mad(vbyu2, mad(vbyu2, -0x1.55cd22p-5f, -0x1.26cf76p-2f), -0x1.55554ep-2f) * vbyu2 * vbyu;
  53     float q = mad(vbyu2, mad(vbyu2, mad(vbyu2, 0x1.9f1304p-5f, 0x1.2656fap-1f), 0x1.76b4b8p+0f), 1.0f);
  54 #endif
  55
  56     // Octant 0 result
  57     float a = mad(p, MATH_RECIP(q), vbyu);
  58
  59     // Fix up 3 other octants
  60     float at = piby2 - a;
  61     a = ay > ax ? at : a;
  62     at = pi - a;
  63     a = x < 0.0F ? at : a;
  64
  65     // y == 0 => 0 for x >= 0, pi for x < 0
  66     at = as_int(x) < 0 ? pi : 0.0f;
  67     a = y == 0.0f ? at : a;
  68
  69     // if (!FINITE_ONLY()) {
  70         // x and y are +- Inf
  71         at = x > 0.0f ? piby4 : threepiby4;
  72         a = ax == INFINITY & ay == INFINITY ? at : a;
  73
  74         // x or y is NaN
  75         a = isnan(x) | isnan(y) ? as_float(QNANBITPATT_SP32) : a;
  76     // }
  77
  78     // Fixup sign and return
  79     return copysign(a, y);
  80 }
  81
  82 _CLC_BINARY_VECTORIZE(_CLC_OVERLOAD _CLC_DEF, float, atan2, float, float);
  83
  84 #ifdef cl_khr_fp64
  85
  86 #pragma OPENCL EXTENSION cl_khr_fp64 : enable
  87
  88 _CLC_OVERLOAD _CLC_DEF double atan2(double y, double x)
  89 {
  90     const double pi = 3.1415926535897932e+00;          /* 0x400921fb54442d18 */
  91     const double piby2 = 1.5707963267948966e+00;       /* 0x3ff921fb54442d18 */
  92     const double piby4 = 7.8539816339744831e-01;       /* 0x3fe921fb54442d18 */
  93     const double three_piby4 = 2.3561944901923449e+00; /* 0x4002d97c7f3321d2 */
  94     const double pi_head = 3.1415926218032836e+00;     /* 0x400921fb50000000 */
  95     const double pi_tail = 3.1786509547056392e-08;     /* 0x3e6110b4611a6263 */
  96     const double piby2_head = 1.5707963267948965e+00;  /* 0x3ff921fb54442d18 */
  97     const double piby2_tail = 6.1232339957367660e-17;  /* 0x3c91a62633145c07 */
  98
  99     double x2 = x;
 100     int xneg = as_int2(x).hi < 0;
 101     int xexp = (as_int2(x).hi >> 20) & 0x7ff;
 102
 103     double y2 = y;
 104     int yneg = as_int2(y).hi < 0;
 105     int yexp = (as_int2(y).hi >> 20) & 0x7ff;
 106
 107     int cond2 = (xexp < 1021) & (yexp < 1021);
 108     int diffexp = yexp - xexp;
 109
 110     // Scale up both x and y if they are both below 1/4
 111     double x1 = ldexp(x, 1024);
 112     int xexp1 = (as_int2(x1).hi >> 20) & 0x7ff;
 113     double y1 = ldexp(y, 1024);
 114     int yexp1 = (as_int2(y1).hi >> 20) & 0x7ff;
 115     int diffexp1 = yexp1 - xexp1;
 116
 117     diffexp = cond2 ? diffexp1 : diffexp;
 118     x = cond2 ? x1 : x;
 119     y = cond2 ? y1 : y;
 120
 121     // General case: take absolute values of arguments
 122     double u = fabs(x);
 123     double v = fabs(y);
 124
 125     // Swap u and v if necessary to obtain 0 < v < u. Compute v/u.
 126     int swap_vu = u < v;
 127     double uu = u;
 128     u = swap_vu ? v : u;
 129     v = swap_vu ? uu : v;
 130
 131     double vbyu = v / u;
 132     double q1, q2;
 133
 134     // General values of v/u. Use a look-up table and series expansion.
 135
 136     {
 137         double val = vbyu > 0.0625 ? vbyu : 0.063;
 138         int index = convert_int(fma(256.0, val, 0.5));
 139         double2 tv = USE_TABLE(atan_jby256_tbl, index - 16);
 140         q1 = tv.s0;
 141         q2 = tv.s1;
 142         double c = (double)index * 0x1.0p-8;
 143
 144         // We're going to scale u and v by 2^(-u_exponent) to bring them close to 1
 145         // u_exponent could be EMAX so we have to do it in 2 steps
 146         int m = -((int)(as_ulong(u) >> EXPSHIFTBITS_DP64) - EXPBIAS_DP64);
 147         //double um = __amdil_ldexp_f64(u, m);
 148         //double vm = __amdil_ldexp_f64(v, m);
 149         double um = ldexp(u, m);
 150         double vm = ldexp(v, m);
 151
 152         // 26 leading bits of u
 153         double u1 = as_double(as_ulong(um) & 0xfffffffff8000000UL);
 154         double u2 = um - u1;
 155
 156         double r = MATH_DIVIDE(fma(-c, u2, fma(-c, u1, vm)), fma(c, vm, um));
 157
 158         // Polynomial approximation to atan(r)
 159         double s = r * r;
 160         q2 = q2 + fma((s * fma(-s, 0.19999918038989143496, 0.33333333333224095522)), -r, r);
 161     }
 162
 163
 164     double q3, q4;
 165     {
 166         q3 = 0.0;
 167         q4 = vbyu;
 168     }
 169
 170     double q5, q6;
 171     {
 172         double u1 = as_double(as_ulong(u) & 0xffffffff00000000UL);
 173         double u2 = u - u1;
 174         double vu1 = as_double(as_ulong(vbyu) & 0xffffffff00000000UL);
 175         double vu2 = vbyu - vu1;
 176
 177         q5 = 0.0;
 178         double s = vbyu * vbyu;
 179         q6 = vbyu + fma(-vbyu * s,
 180                         fma(-s,
 181                             fma(-s,
 182                                 fma(-s,
 183                                     fma(-s, 0.90029810285449784439E-01,
 184                                         0.11110736283514525407),
 185                                     0.14285713561807169030),
 186                                 0.19999999999393223405),
 187                             0.33333333333333170500),
 188                          MATH_DIVIDE(fma(-u, vu2, fma(-u2, vu1, fma(-u1, vu1, v))), u));
 189     }
 190
 191
 192     q3 = vbyu < 0x1.d12ed0af1a27fp-27 ? q3 : q5;
 193     q4 = vbyu < 0x1.d12ed0af1a27fp-27 ? q4 : q6;
 194
 195     q1 = vbyu > 0.0625 ? q1 : q3;
 196     q2 = vbyu > 0.0625 ? q2 : q4;
 197
 198     // Tidy-up according to which quadrant the arguments lie in
 199     double res1, res2, res3, res4;
 200     q1 = swap_vu ? piby2_head - q1 : q1;
 201     q2 = swap_vu ? piby2_tail - q2 : q2;
 202     q1 = xneg ? pi_head - q1 : q1;
 203     q2 = xneg ? pi_tail - q2 : q2;
 204     q1 = q1 + q2;
 205     res4 = yneg ? -q1 : q1;
 206
 207     res1 = yneg ? -three_piby4 : three_piby4;
 208     res2 = yneg ? -piby4 : piby4;
 209     res3 = xneg ? res1 : res2;
 210
 211     res3 = isinf(x2) & isinf(y2) ? res3 : res4;
 212     res1 = yneg ? -pi : pi;
 213
 214     // abs(x)/abs(y) > 2^56 and x < 0
 215     res3 = (diffexp < -56 && xneg) ? res1 : res3;
 216
 217     res4 = MATH_DIVIDE(y, x);
 218     // x positive and dominant over y by a factor of 2^28
 219     res3 = diffexp < -28 & xneg == 0 ? res4 : res3;
 220
 221     // abs(y)/abs(x) > 2^56
 222     res4 = yneg ? -piby2 : piby2;       // atan(y/x) is insignificant compared to piby2
 223     res3 = diffexp > 56 ? res4 : res3;
 224
 225     res3 = x2 == 0.0 ? res4 : res3;   // Zero x gives +- pi/2 depending on sign of y
 226     res4 = xneg ? res1 : y2;
 227
 228     res3 = y2 == 0.0 ? res4 : res3;   // Zero y gives +-0 for positive x and +-pi for negative x
 229     res3 = isnan(y2) ? y2 : res3;
 230     res3 = isnan(x2) ? x2 : res3;
 231
 232     return res3;
 233 }
 234
 235 _CLC_BINARY_VECTORIZE(_CLC_OVERLOAD _CLC_DEF, double, atan2, double, double);
 236
 237 #endif