249 lines · plain
1//===----------------------------------------------------------------------===//2//3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.4// See https://llvm.org/LICENSE.txt for license information.5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception6//7//===----------------------------------------------------------------------===//8 9#if __CLC_FPSIZE == 3210 11_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_atan2(__CLC_GENTYPE y,12 __CLC_GENTYPE x) {13 const __CLC_GENTYPE pi = 0x1.921fb6p+1f;14 const __CLC_GENTYPE piby2 = 0x1.921fb6p+0f;15 const __CLC_GENTYPE piby4 = 0x1.921fb6p-1f;16 const __CLC_GENTYPE threepiby4 = 0x1.2d97c8p+1f;17 18 __CLC_GENTYPE ax = __clc_fabs(x);19 __CLC_GENTYPE ay = __clc_fabs(y);20 __CLC_GENTYPE v = __clc_min(ax, ay);21 __CLC_GENTYPE u = __clc_max(ax, ay);22 23 // Scale since u could be large, as in "regular" divide24 __CLC_GENTYPE s = u > 0x1.0p+96f ? 0x1.0p-32f : 1.0f;25 __CLC_GENTYPE vbyu = s * MATH_DIVIDE(v, s * u);26 27 __CLC_GENTYPE vbyu2 = vbyu * vbyu;28 29#define USE_2_2_APPROXIMATION30#if defined USE_2_2_APPROXIMATION31 __CLC_GENTYPE p =32 __clc_mad(vbyu2, __clc_mad(vbyu2, -0x1.7e1f78p-9f, -0x1.7d1b98p-3f),33 -0x1.5554d0p-2f) *34 vbyu2 * vbyu;35 __CLC_GENTYPE q =36 __clc_mad(vbyu2, __clc_mad(vbyu2, 0x1.1a714cp-2f, 0x1.287c56p+0f), 1.0f);37#else38 __CLC_GENTYPE p =39 __clc_mad(vbyu2, __clc_mad(vbyu2, -0x1.55cd22p-5f, -0x1.26cf76p-2f),40 -0x1.55554ep-2f) *41 vbyu2 * vbyu;42 __CLC_GENTYPE q = __clc_mad(43 vbyu2,44 __clc_mad(vbyu2, __clc_mad(vbyu2, 0x1.9f1304p-5f, 0x1.2656fap-1f),45 0x1.76b4b8p+0f),46 1.0f);47#endif48 49 // Octant 0 result50 __CLC_GENTYPE a = __clc_mad(p, MATH_RECIP(q), vbyu);51 52 // Fix up 3 other octants53 __CLC_GENTYPE at = piby2 - a;54 a = ay > ax ? at : a;55 at = pi - a;56 a = x < 0.0F ? at : a;57 58 // y == 0 => 0 for x >= 0, pi for x < 059 at = __CLC_AS_INTN(x) < 0 ? pi : 0.0f;60 a = y == 0.0f ? at : a;61 62 // x and y are +- Inf63 at = x > 0.0f ? piby4 : threepiby4;64 a = __clc_select(a, at, __clc_isinf(x) && __clc_isinf(y));65 66 // x or y is NaN67 a = __clc_select(a, __CLC_GENTYPE_NAN, __clc_isnan(x) || __clc_isnan(y));68 69 // Fixup sign and return70 return __clc_copysign(a, y);71}72 73#elif __CLC_FPSIZE == 6474 75_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_atan2(__CLC_GENTYPE y,76 __CLC_GENTYPE x) {77 const __CLC_GENTYPE pi = 3.1415926535897932e+00; /* 0x400921fb54442d18 */78 const __CLC_GENTYPE piby2 = 1.5707963267948966e+00; /* 0x3ff921fb54442d18 */79 const __CLC_GENTYPE piby4 = 7.8539816339744831e-01; /* 0x3fe921fb54442d18 */80 // 0x4002d97c7f3321d281 const __CLC_GENTYPE three_piby4 = 2.3561944901923449e+00;82 const __CLC_GENTYPE pi_head = 3.1415926218032836e+00; /* 0x400921fb50000000 */83 const __CLC_GENTYPE pi_tail = 3.1786509547056392e-08; /* 0x3e6110b4611a6263 */84 // 0x3ff921fb54442d1885 const __CLC_GENTYPE piby2_head = 1.5707963267948965e+00;86 // 0x3c91a62633145c0787 const __CLC_GENTYPE piby2_tail = 6.1232339957367660e-17;88 89 __CLC_GENTYPE x2 = x;90 // Important to capture -0.0 in xneg and yneg, so comparison done as integer91 __CLC_LONGN xneg = __CLC_AS_LONGN(x) < 0;92 __CLC_INTN xexp =93 __CLC_CONVERT_INTN(__CLC_AS_ULONGN(x) >> EXPSHIFTBITS_DP64) & 0x7ff;94 95 __CLC_GENTYPE y2 = y;96 __CLC_LONGN yneg = __CLC_AS_LONGN(y) < 0;97 __CLC_INTN yexp =98 __CLC_CONVERT_INTN(__CLC_AS_ULONGN(y) >> EXPSHIFTBITS_DP64) & 0x7ff;99 100 __CLC_LONGN cond2 = __CLC_CONVERT_LONGN(xexp < 1021 && yexp < 1021);101 __CLC_LONGN diffexp = __CLC_CONVERT_LONGN(yexp - xexp);102 103 // Scale up both x and y if they are both below 1/4104 __CLC_GENTYPE x1 = __clc_ldexp(x, 1024);105 __CLC_INTN xexp1 =106 __CLC_CONVERT_INTN(__CLC_AS_ULONGN(x1) >> EXPSHIFTBITS_DP64) & 0x7ff;107 __CLC_GENTYPE y1 = __clc_ldexp(y, 1024);108 __CLC_INTN yexp1 =109 __CLC_CONVERT_INTN(__CLC_AS_ULONGN(y1) >> EXPSHIFTBITS_DP64) & 0x7ff;110 __CLC_LONGN diffexp1 = __CLC_CONVERT_LONGN(yexp1 - xexp1);111 112 diffexp = __clc_select(diffexp, diffexp1, cond2);113 x = cond2 ? x1 : x;114 y = cond2 ? y1 : y;115 116 // General case: take absolute values of arguments117 __CLC_GENTYPE u = __clc_fabs(x);118 __CLC_GENTYPE v = __clc_fabs(y);119 120 // Swap u and v if necessary to obtain 0 < v < u. Compute v/u.121 __CLC_LONGN swap_vu = u < v;122 __CLC_GENTYPE uu = u;123 u = swap_vu ? v : u;124 v = swap_vu ? uu : v;125 126 __CLC_GENTYPE vbyu = v / u;127 __CLC_GENTYPE q1, q2;128 129 // General values of v/u. Use a look-up table and series expansion.130 131 {132 __CLC_GENTYPE val = vbyu > 0.0625 ? vbyu : 0.063;133 __CLC_INTN index = __CLC_CONVERT_INTN(__clc_fma(256.0, val, 0.5));134 q1 = __CLC_USE_TABLE(atan_jby256_tbl_head, index - 16);135 q2 = __CLC_USE_TABLE(atan_jby256_tbl_tail, index - 16);136 __CLC_GENTYPE c = __CLC_CONVERT_GENTYPE(index) * 0x1.0p-8;137 138 // We're going to scale u and v by 2^(-u_exponent) to bring them close to 1139 // u_exponent could be EMAX so we have to do it in 2 steps140 __CLC_INTN m =141 -(__CLC_CONVERT_INTN(__CLC_AS_ULONGN(u) >> EXPSHIFTBITS_DP64) -142 EXPBIAS_DP64);143 __CLC_GENTYPE um = __clc_ldexp(u, m);144 __CLC_GENTYPE vm = __clc_ldexp(v, m);145 146 // 26 leading bits of u147 __CLC_GENTYPE u1 =148 __CLC_AS_GENTYPE(__CLC_AS_ULONGN(um) & 0xfffffffff8000000UL);149 __CLC_GENTYPE u2 = um - u1;150 151 __CLC_GENTYPE r = MATH_DIVIDE(__clc_fma(-c, u2, __clc_fma(-c, u1, vm)),152 __clc_fma(c, vm, um));153 154 // Polynomial approximation to atan(r)155 __CLC_GENTYPE s = r * r;156 q2 = q2 + __clc_fma((s * __clc_fma(-s, 0.19999918038989143496,157 0.33333333333224095522)),158 -r, r);159 }160 161 __CLC_GENTYPE q3, q4;162 {163 q3 = 0.0;164 q4 = vbyu;165 }166 167 __CLC_GENTYPE q5, q6;168 {169 __CLC_GENTYPE u1 =170 __CLC_AS_GENTYPE(__CLC_AS_ULONGN(u) & 0xffffffff00000000UL);171 __CLC_GENTYPE u2 = u - u1;172 __CLC_GENTYPE vu1 =173 __CLC_AS_GENTYPE(__CLC_AS_ULONGN(vbyu) & 0xffffffff00000000UL);174 __CLC_GENTYPE vu2 = vbyu - vu1;175 176 q5 = 0.0;177 __CLC_GENTYPE s = vbyu * vbyu;178 q6 = vbyu +179 __clc_fma(180 -vbyu * s,181 __clc_fma(182 -s,183 __clc_fma(-s,184 __clc_fma(-s,185 __clc_fma(-s, 0.90029810285449784439E-01,186 0.11110736283514525407),187 0.14285713561807169030),188 0.19999999999393223405),189 0.33333333333333170500),190 MATH_DIVIDE(__clc_fma(-u, vu2,191 __clc_fma(-u2, vu1, __clc_fma(-u1, vu1, v))),192 u));193 }194 195 q3 = vbyu < 0x1.d12ed0af1a27fp-27 ? q3 : q5;196 q4 = vbyu < 0x1.d12ed0af1a27fp-27 ? q4 : q6;197 198 q1 = vbyu > 0.0625 ? q1 : q3;199 q2 = vbyu > 0.0625 ? q2 : q4;200 201 // Tidy-up according to which quadrant the arguments lie in202 __CLC_GENTYPE res1, res2, res3, res4;203 q1 = swap_vu ? piby2_head - q1 : q1;204 q2 = swap_vu ? piby2_tail - q2 : q2;205 q1 = xneg ? pi_head - q1 : q1;206 q2 = xneg ? pi_tail - q2 : q2;207 q1 = q1 + q2;208 res4 = yneg ? -q1 : q1;209 210 res1 = yneg ? -three_piby4 : three_piby4;211 res2 = yneg ? -piby4 : piby4;212 res3 = xneg ? res1 : res2;213 214 res3 = __clc_select(res4, res3,215 __CLC_CONVERT_LONGN(__clc_isinf(x2) && __clc_isinf(y2)));216 res1 = yneg ? -pi : pi;217 218 // abs(x)/abs(y) > 2^56 and x < 0219 res3 = (diffexp < -56 && xneg) ? res1 : res3;220 221 res4 = MATH_DIVIDE(y, x);222 // x positive and dominant over y by a factor of 2^28223 res3 = diffexp < -28 && xneg == 0 ? res4 : res3;224 225 // abs(y)/abs(x) > 2^56226 res4 = yneg ? -piby2 : piby2; // atan(y/x) is insignificant compared to piby2227 res3 = diffexp > 56 ? res4 : res3;228 229 res3 = x2 == 0.0 ? res4 : res3; // Zero x gives +- pi/2 depending on sign of y230 res4 = xneg ? res1 : y2;231 232 // Zero y gives +-0 for positive x and +-pi for negative x233 res3 = y2 == 0.0 ? res4 : res3;234 res3 = __clc_isnan(y2) ? y2 : res3;235 res3 = __clc_isnan(x2) ? x2 : res3;236 237 return res3;238}239 240#elif __CLC_FPSIZE == 16241 242_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_atan2(__CLC_GENTYPE x,243 __CLC_GENTYPE y) {244 return __CLC_CONVERT_GENTYPE(245 __clc_atan2(__CLC_CONVERT_FLOATN(x), __CLC_CONVERT_FLOATN(y)));246}247 248#endif249