228 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_atan2pi(__CLC_GENTYPE y,12 __CLC_GENTYPE x) {13 const __CLC_GENTYPE pi = 0x1.921fb6p+1f;14 15 __CLC_GENTYPE ax = __clc_fabs(x);16 __CLC_GENTYPE ay = __clc_fabs(y);17 __CLC_GENTYPE v = __clc_min(ax, ay);18 __CLC_GENTYPE u = __clc_max(ax, ay);19 20 // Scale since u could be large, as in "regular" divide21 __CLC_GENTYPE s = u > 0x1.0p+96f ? 0x1.0p-32f : 1.0f;22 __CLC_GENTYPE vbyu = s * MATH_DIVIDE(v, s * u);23 24 __CLC_GENTYPE vbyu2 = vbyu * vbyu;25 26 __CLC_GENTYPE p =27 __clc_mad(vbyu2, __clc_mad(vbyu2, -0x1.7e1f78p-9f, -0x1.7d1b98p-3f),28 -0x1.5554d0p-2f) *29 vbyu2 * vbyu;30 __CLC_GENTYPE q =31 __clc_mad(vbyu2, __clc_mad(vbyu2, 0x1.1a714cp-2f, 0x1.287c56p+0f), 1.0f);32 33 // Octant 0 result34 __CLC_GENTYPE a = MATH_DIVIDE(__clc_mad(p, MATH_RECIP(q), vbyu), pi);35 36 // Fix up 3 other octants37 __CLC_GENTYPE at = 0.5f - a;38 a = ay > ax ? at : a;39 at = 1.0f - a;40 a = x < 0.0F ? at : a;41 42 // y == 0 => 0 for x >= 0, pi for x < 043 at = __CLC_AS_INTN(x) < 0 ? 1.0f : 0.0f;44 a = y == 0.0f ? at : a;45 46 // x and y are +- Inf47 at = x > 0.0f ? 0.25f : 0.75f;48 a = __clc_select(a, at, __clc_isinf(x) && __clc_isinf(y));49 50 // x or y is NaN51 a = __clc_select(a, __CLC_GENTYPE_NAN, __clc_isnan(x) || __clc_isnan(y));52 53 // Fixup sign and return54 return __clc_copysign(a, y);55}56 57#elif __CLC_FPSIZE == 6458 59_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_atan2pi(__CLC_GENTYPE y,60 __CLC_GENTYPE x) {61 const __CLC_GENTYPE pi = 3.1415926535897932e+00; /* 0x400921fb54442d18 */62 const __CLC_GENTYPE pi_head = 3.1415926218032836e+00; /* 0x400921fb50000000 */63 const __CLC_GENTYPE pi_tail = 3.1786509547056392e-08; /* 0x3e6110b4611a6263 */64 // 0x3ff921fb54442d1865 const __CLC_GENTYPE piby2_head = 1.5707963267948965e+00;66 // 0x3c91a62633145c0767 const __CLC_GENTYPE piby2_tail = 6.1232339957367660e-17;68 69 __CLC_GENTYPE x2 = x;70 __CLC_LONGN xneg = __CLC_AS_LONGN(x) < 0;71 __CLC_INTN xexp =72 __CLC_CONVERT_INTN(__CLC_AS_ULONGN(x) >> EXPSHIFTBITS_DP64) & 0x7ff;73 74 __CLC_GENTYPE y2 = y;75 __CLC_LONGN yneg = __CLC_AS_LONGN(y) < 0;76 __CLC_INTN yexp =77 __CLC_CONVERT_INTN(__CLC_AS_ULONGN(y) >> EXPSHIFTBITS_DP64) & 0x7ff;78 79 __CLC_LONGN cond2 = __CLC_CONVERT_LONGN(xexp < 1021 & yexp < 1021);80 __CLC_LONGN diffexp = __CLC_CONVERT_LONGN(yexp - xexp);81 82 // Scale up both x and y if they are both below 1/483 __CLC_GENTYPE x1 = __clc_ldexp(x, 1024);84 __CLC_INTN xexp1 =85 __CLC_CONVERT_INTN(__CLC_AS_ULONGN(x1) >> EXPSHIFTBITS_DP64) & 0x7ff;86 __CLC_GENTYPE y1 = __clc_ldexp(y, 1024);87 __CLC_INTN yexp1 =88 __CLC_CONVERT_INTN(__CLC_AS_ULONGN(y1) >> EXPSHIFTBITS_DP64) & 0x7ff;89 __CLC_LONGN diffexp1 = __CLC_CONVERT_LONGN(yexp1 - xexp1);90 91 diffexp = __clc_select(diffexp, diffexp1, cond2);92 x = cond2 ? x1 : x;93 y = cond2 ? y1 : y;94 95 // General case: take absolute values of arguments96 __CLC_GENTYPE u = __clc_fabs(x);97 __CLC_GENTYPE v = __clc_fabs(y);98 99 // Swap u and v if necessary to obtain 0 < v < u. Compute v/u.100 __CLC_LONGN swap_vu = u < v;101 __CLC_GENTYPE uu = u;102 u = swap_vu ? v : u;103 v = swap_vu ? uu : v;104 105 __CLC_GENTYPE vbyu = v / u;106 __CLC_GENTYPE q1, q2;107 108 // General values of v/u. Use a look-up table and series expansion.109 110 {111 __CLC_GENTYPE val = vbyu > 0.0625 ? vbyu : 0.063;112 __CLC_INTN index = __CLC_CONVERT_INTN(__clc_fma(256.0, val, 0.5));113 q1 = __CLC_USE_TABLE(atan_jby256_tbl_head, (index - 16));114 q2 = __CLC_USE_TABLE(atan_jby256_tbl_tail, (index - 16));115 __CLC_GENTYPE c = __CLC_CONVERT_GENTYPE(index) * 0x1.0p-8;116 117 // We're going to scale u and v by 2^(-u_exponent) to bring them close to 1118 // u_exponent could be EMAX so we have to do it in 2 steps119 __CLC_INTN m =120 -(__CLC_CONVERT_INTN(__CLC_AS_ULONGN(u) >> EXPSHIFTBITS_DP64) -121 EXPBIAS_DP64);122 __CLC_GENTYPE um = __clc_ldexp(u, m);123 __CLC_GENTYPE vm = __clc_ldexp(v, m);124 125 // 26 leading bits of u126 __CLC_GENTYPE u1 =127 __CLC_AS_GENTYPE(__CLC_AS_ULONGN(um) & 0xfffffffff8000000UL);128 __CLC_GENTYPE u2 = um - u1;129 130 __CLC_GENTYPE r = MATH_DIVIDE(__clc_fma(-c, u2, __clc_fma(-c, u1, vm)),131 __clc_fma(c, vm, um));132 133 // Polynomial approximation to atan(r)134 __CLC_GENTYPE s = r * r;135 q2 = q2 + __clc_fma((s * __clc_fma(-s, 0.19999918038989143496,136 0.33333333333224095522)),137 -r, r);138 }139 140 __CLC_GENTYPE q3, q4;141 {142 q3 = 0.0;143 q4 = vbyu;144 }145 146 __CLC_GENTYPE q5, q6;147 {148 __CLC_GENTYPE u1 =149 __CLC_AS_GENTYPE(__CLC_AS_ULONGN(u) & 0xffffffff00000000UL);150 __CLC_GENTYPE u2 = u - u1;151 __CLC_GENTYPE vu1 =152 __CLC_AS_GENTYPE(__CLC_AS_ULONGN(vbyu) & 0xffffffff00000000UL);153 __CLC_GENTYPE vu2 = vbyu - vu1;154 155 q5 = 0.0;156 __CLC_GENTYPE s = vbyu * vbyu;157 q6 = vbyu +158 __clc_fma(159 -vbyu * s,160 __clc_fma(161 -s,162 __clc_fma(-s,163 __clc_fma(-s,164 __clc_fma(-s, 0.90029810285449784439E-01,165 0.11110736283514525407),166 0.14285713561807169030),167 0.19999999999393223405),168 0.33333333333333170500),169 MATH_DIVIDE(__clc_fma(-u, vu2,170 __clc_fma(-u2, vu1, __clc_fma(-u1, vu1, v))),171 u));172 }173 174 q3 = vbyu < 0x1.d12ed0af1a27fp-27 ? q3 : q5;175 q4 = vbyu < 0x1.d12ed0af1a27fp-27 ? q4 : q6;176 177 q1 = vbyu > 0.0625 ? q1 : q3;178 q2 = vbyu > 0.0625 ? q2 : q4;179 180 // Tidy-up according to which quadrant the arguments lie in181 __CLC_GENTYPE res1, res2, res3, res4;182 q1 = swap_vu ? piby2_head - q1 : q1;183 q2 = swap_vu ? piby2_tail - q2 : q2;184 q1 = xneg ? pi_head - q1 : q1;185 q2 = xneg ? pi_tail - q2 : q2;186 q1 = MATH_DIVIDE(q1 + q2, pi);187 res4 = yneg ? -q1 : q1;188 189 res1 = yneg ? -0.75 : 0.75;190 res2 = yneg ? -0.25 : 0.25;191 res3 = xneg ? res1 : res2;192 193 res3 = __clc_select(res4, res3,194 __CLC_CONVERT_LONGN(__clc_isinf(y2) & __clc_isinf(x2)));195 res1 = yneg ? -1.0 : 1.0;196 197 // abs(x)/abs(y) > 2^56 and x < 0198 res3 = diffexp < -56 && xneg ? res1 : res3;199 200 res4 = MATH_DIVIDE(MATH_DIVIDE(y, x), pi);201 // x positive and dominant over y by a factor of 2^28202 res3 = diffexp < -28 && xneg == 0 ? res4 : res3;203 204 // abs(y)/abs(x) > 2^56205 res4 = yneg ? -0.5 : 0.5; // atan(y/x) is insignificant compared to piby2206 res3 = diffexp > 56 ? res4 : res3;207 208 res3 = x2 == 0.0 ? res4 : res3; // Zero x gives +- pi/2 depending on sign of y209 res4 = xneg ? res1 : y2;210 211 // Zero y gives +-0 for positive x and +-pi for negative x212 res3 = y2 == 0.0 ? res4 : res3;213 res3 = __clc_isnan(y2) ? y2 : res3;214 res3 = __clc_isnan(x2) ? x2 : res3;215 216 return res3;217}218 219#elif __CLC_FPSIZE == 16220 221_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_atan2pi(__CLC_GENTYPE x,222 __CLC_GENTYPE y) {223 return __CLC_CONVERT_GENTYPE(224 __clc_atan2pi(__CLC_CONVERT_FLOATN(x), __CLC_CONVERT_FLOATN(y)));225}226 227#endif228