libclc: Update atan2pi (#188707)
This was originally ported from rocm device libs in 37406a209c75a09f850cd5e5498568d34a6f05d1. Merge in more recent changes.
This commit is contained in:
Matt Arsenault
GitHub
parent
207598a827
commit
8430e9e8f0
@ -6,19 +6,19 @@
|
||||
//
|
||||
//===----------------------------------------------------------------------===//
|
||||
|
||||
#include "clc/clc_convert.h"
|
||||
#include "clc/float/definitions.h"
|
||||
#include "clc/internal/clc.h"
|
||||
#include "clc/math/clc_atan_helpers.h"
|
||||
#include "clc/math/clc_copysign.h"
|
||||
#include "clc/math/clc_fabs.h"
|
||||
#include "clc/math/clc_fma.h"
|
||||
#include "clc/math/clc_ldexp.h"
|
||||
#include "clc/math/clc_fmax.h"
|
||||
#include "clc/math/clc_fmin.h"
|
||||
#include "clc/math/clc_mad.h"
|
||||
#include "clc/math/math.h"
|
||||
#include "clc/math/tables.h"
|
||||
#include "clc/relational/clc_isinf.h"
|
||||
#include "clc/relational/clc_isnan.h"
|
||||
#include "clc/relational/clc_isunordered.h"
|
||||
#include "clc/relational/clc_select.h"
|
||||
#include "clc/relational/clc_signbit.h"
|
||||
#include "clc/shared/clc_max.h"
|
||||
#include "clc/shared/clc_min.h"
|
||||
|
||||
|
||||
@ -6,222 +6,99 @@
|
||||
//
|
||||
//===----------------------------------------------------------------------===//
|
||||
|
||||
#pragma OPENCL FP_CONTRACT OFF
|
||||
|
||||
#if __CLC_FPSIZE == 32
|
||||
|
||||
_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_atan2pi(__CLC_GENTYPE y,
|
||||
__CLC_GENTYPE x) {
|
||||
const __CLC_GENTYPE pi = 0x1.921fb6p+1f;
|
||||
_CLC_DEF _CLC_OVERLOAD _CLC_CONST __CLC_FLOATN __clc_atan2pi(__CLC_FLOATN y,
|
||||
__CLC_FLOATN x) {
|
||||
__CLC_FLOATN ax = __clc_fabs(x);
|
||||
__CLC_FLOATN ay = __clc_fabs(y);
|
||||
__CLC_FLOATN v = __clc_fmin(ax, ay);
|
||||
__CLC_FLOATN u = __clc_fmax(ax, ay);
|
||||
|
||||
__CLC_GENTYPE ax = __clc_fabs(x);
|
||||
__CLC_GENTYPE ay = __clc_fabs(y);
|
||||
__CLC_GENTYPE v = __clc_min(ax, ay);
|
||||
__CLC_GENTYPE u = __clc_max(ax, ay);
|
||||
__CLC_FLOATN vbyu = v / u;
|
||||
|
||||
// Scale since u could be large, as in "regular" divide
|
||||
__CLC_GENTYPE s = u > 0x1.0p+96f ? 0x1.0p-32f : 1.0f;
|
||||
__CLC_GENTYPE vbyu = s * MATH_DIVIDE(v, s * u);
|
||||
__CLC_FLOATN a = __clc_atanpi_reduced(vbyu);
|
||||
|
||||
__CLC_GENTYPE vbyu2 = vbyu * vbyu;
|
||||
|
||||
__CLC_GENTYPE p =
|
||||
__clc_mad(vbyu2, __clc_mad(vbyu2, -0x1.7e1f78p-9f, -0x1.7d1b98p-3f),
|
||||
-0x1.5554d0p-2f) *
|
||||
vbyu2 * vbyu;
|
||||
__CLC_GENTYPE q =
|
||||
__clc_mad(vbyu2, __clc_mad(vbyu2, 0x1.1a714cp-2f, 0x1.287c56p+0f), 1.0f);
|
||||
|
||||
// Octant 0 result
|
||||
__CLC_GENTYPE a = MATH_DIVIDE(__clc_mad(p, MATH_RECIP(q), vbyu), pi);
|
||||
|
||||
// Fix up 3 other octants
|
||||
__CLC_GENTYPE at = 0.5f - a;
|
||||
__CLC_FLOATN at = 0.5f - a;
|
||||
a = ay > ax ? at : a;
|
||||
at = 1.0f - a;
|
||||
a = x < 0.0F ? at : a;
|
||||
a = x < 0.0f ? at : a;
|
||||
|
||||
// y == 0 => 0 for x >= 0, pi for x < 0
|
||||
at = __CLC_AS_INTN(x) < 0 ? 1.0f : 0.0f;
|
||||
at = __clc_signbit(x) ? 1.0f : 0.0f;
|
||||
a = y == 0.0f ? at : a;
|
||||
|
||||
// x and y are +- Inf
|
||||
at = x > 0.0f ? 0.25f : 0.75f;
|
||||
a = __clc_select(a, at, __clc_isinf(x) && __clc_isinf(y));
|
||||
at = x < 0.0f ? 0.75f : 0.25f;
|
||||
a = (__clc_isinf(x) & __clc_isinf(y)) ? at : a;
|
||||
|
||||
// x or y is NaN
|
||||
a = __clc_select(a, __CLC_GENTYPE_NAN, __clc_isnan(x) || __clc_isnan(y));
|
||||
a = __clc_isunordered(x, y) ? FLT_NAN : a;
|
||||
|
||||
// Fixup sign and return
|
||||
return __clc_copysign(a, y);
|
||||
}
|
||||
|
||||
#elif __CLC_FPSIZE == 64
|
||||
|
||||
_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_atan2pi(__CLC_GENTYPE y,
|
||||
__CLC_GENTYPE x) {
|
||||
const __CLC_GENTYPE pi = 3.1415926535897932e+00; /* 0x400921fb54442d18 */
|
||||
const __CLC_GENTYPE pi_head = 3.1415926218032836e+00; /* 0x400921fb50000000 */
|
||||
const __CLC_GENTYPE pi_tail = 3.1786509547056392e-08; /* 0x3e6110b4611a6263 */
|
||||
// 0x3ff921fb54442d18
|
||||
const __CLC_GENTYPE piby2_head = 1.5707963267948965e+00;
|
||||
// 0x3c91a62633145c07
|
||||
const __CLC_GENTYPE piby2_tail = 6.1232339957367660e-17;
|
||||
_CLC_DEF _CLC_OVERLOAD _CLC_CONST __CLC_DOUBLEN __clc_atan2pi(__CLC_DOUBLEN y,
|
||||
__CLC_DOUBLEN x) {
|
||||
__CLC_DOUBLEN ay = __clc_fabs(y);
|
||||
__CLC_DOUBLEN ax = __clc_fabs(x);
|
||||
__CLC_DOUBLEN u = __clc_fmax(ax, ay);
|
||||
__CLC_DOUBLEN v = __clc_fmin(ax, ay);
|
||||
__CLC_DOUBLEN vbyu = v / u;
|
||||
|
||||
__CLC_GENTYPE x2 = x;
|
||||
__CLC_LONGN xneg = __CLC_AS_LONGN(x) < 0;
|
||||
__CLC_INTN xexp =
|
||||
__CLC_CONVERT_INTN(__CLC_AS_ULONGN(x) >> EXPSHIFTBITS_DP64) & 0x7ff;
|
||||
__CLC_DOUBLEN a = __clc_atanpi_reduced(vbyu);
|
||||
|
||||
__CLC_GENTYPE y2 = y;
|
||||
__CLC_LONGN yneg = __CLC_AS_LONGN(y) < 0;
|
||||
__CLC_INTN yexp =
|
||||
__CLC_CONVERT_INTN(__CLC_AS_ULONGN(y) >> EXPSHIFTBITS_DP64) & 0x7ff;
|
||||
__CLC_LONGN xneg = __clc_signbit(x);
|
||||
|
||||
__CLC_LONGN cond2 = __CLC_CONVERT_LONGN(xexp < 1021 & yexp < 1021);
|
||||
__CLC_LONGN diffexp = __CLC_CONVERT_LONGN(yexp - xexp);
|
||||
__CLC_DOUBLEN t = 0.5 - a;
|
||||
a = ax < ay ? t : a;
|
||||
t = 1.0 - a;
|
||||
a = xneg ? t : a;
|
||||
|
||||
// Scale up both x and y if they are both below 1/4
|
||||
__CLC_GENTYPE x1 = __clc_ldexp(x, 1024);
|
||||
__CLC_INTN xexp1 =
|
||||
__CLC_CONVERT_INTN(__CLC_AS_ULONGN(x1) >> EXPSHIFTBITS_DP64) & 0x7ff;
|
||||
__CLC_GENTYPE y1 = __clc_ldexp(y, 1024);
|
||||
__CLC_INTN yexp1 =
|
||||
__CLC_CONVERT_INTN(__CLC_AS_ULONGN(y1) >> EXPSHIFTBITS_DP64) & 0x7ff;
|
||||
__CLC_LONGN diffexp1 = __CLC_CONVERT_LONGN(yexp1 - xexp1);
|
||||
t = xneg ? 1.0 : 0.0;
|
||||
a = y == 0.0 ? t : a;
|
||||
|
||||
diffexp = __clc_select(diffexp, diffexp1, cond2);
|
||||
x = cond2 ? x1 : x;
|
||||
y = cond2 ? y1 : y;
|
||||
t = xneg ? 0.75 : 0.25;
|
||||
t = __clc_copysign(t, y);
|
||||
a = (__clc_isinf(x) & __clc_isinf(y)) ? t : a;
|
||||
|
||||
// General case: take absolute values of arguments
|
||||
__CLC_GENTYPE u = __clc_fabs(x);
|
||||
__CLC_GENTYPE v = __clc_fabs(y);
|
||||
a = __clc_isunordered(x, y) ? DBL_NAN : a;
|
||||
|
||||
// Swap u and v if necessary to obtain 0 < v < u. Compute v/u.
|
||||
__CLC_LONGN swap_vu = u < v;
|
||||
__CLC_GENTYPE uu = u;
|
||||
u = swap_vu ? v : u;
|
||||
v = swap_vu ? uu : v;
|
||||
|
||||
__CLC_GENTYPE vbyu = v / u;
|
||||
__CLC_GENTYPE q1, q2;
|
||||
|
||||
// General values of v/u. Use a look-up table and series expansion.
|
||||
|
||||
{
|
||||
__CLC_GENTYPE val = vbyu > 0.0625 ? vbyu : 0.063;
|
||||
__CLC_INTN index = __CLC_CONVERT_INTN(__clc_fma(256.0, val, 0.5));
|
||||
q1 = __CLC_USE_TABLE(atan_jby256_tbl_head, (index - 16));
|
||||
q2 = __CLC_USE_TABLE(atan_jby256_tbl_tail, (index - 16));
|
||||
__CLC_GENTYPE c = __CLC_CONVERT_GENTYPE(index) * 0x1.0p-8;
|
||||
|
||||
// We're going to scale u and v by 2^(-u_exponent) to bring them close to 1
|
||||
// u_exponent could be EMAX so we have to do it in 2 steps
|
||||
__CLC_INTN m =
|
||||
-(__CLC_CONVERT_INTN(__CLC_AS_ULONGN(u) >> EXPSHIFTBITS_DP64) -
|
||||
EXPBIAS_DP64);
|
||||
__CLC_GENTYPE um = __clc_ldexp(u, m);
|
||||
__CLC_GENTYPE vm = __clc_ldexp(v, m);
|
||||
|
||||
// 26 leading bits of u
|
||||
__CLC_GENTYPE u1 =
|
||||
__CLC_AS_GENTYPE(__CLC_AS_ULONGN(um) & 0xfffffffff8000000UL);
|
||||
__CLC_GENTYPE u2 = um - u1;
|
||||
|
||||
__CLC_GENTYPE r = MATH_DIVIDE(__clc_fma(-c, u2, __clc_fma(-c, u1, vm)),
|
||||
__clc_fma(c, vm, um));
|
||||
|
||||
// Polynomial approximation to atan(r)
|
||||
__CLC_GENTYPE s = r * r;
|
||||
q2 = q2 + __clc_fma((s * __clc_fma(-s, 0.19999918038989143496,
|
||||
0.33333333333224095522)),
|
||||
-r, r);
|
||||
}
|
||||
|
||||
__CLC_GENTYPE q3, q4;
|
||||
{
|
||||
q3 = 0.0;
|
||||
q4 = vbyu;
|
||||
}
|
||||
|
||||
__CLC_GENTYPE q5, q6;
|
||||
{
|
||||
__CLC_GENTYPE u1 =
|
||||
__CLC_AS_GENTYPE(__CLC_AS_ULONGN(u) & 0xffffffff00000000UL);
|
||||
__CLC_GENTYPE u2 = u - u1;
|
||||
__CLC_GENTYPE vu1 =
|
||||
__CLC_AS_GENTYPE(__CLC_AS_ULONGN(vbyu) & 0xffffffff00000000UL);
|
||||
__CLC_GENTYPE vu2 = vbyu - vu1;
|
||||
|
||||
q5 = 0.0;
|
||||
__CLC_GENTYPE s = vbyu * vbyu;
|
||||
q6 = vbyu +
|
||||
__clc_fma(
|
||||
-vbyu * s,
|
||||
__clc_fma(
|
||||
-s,
|
||||
__clc_fma(-s,
|
||||
__clc_fma(-s,
|
||||
__clc_fma(-s, 0.90029810285449784439E-01,
|
||||
0.11110736283514525407),
|
||||
0.14285713561807169030),
|
||||
0.19999999999393223405),
|
||||
0.33333333333333170500),
|
||||
MATH_DIVIDE(__clc_fma(-u, vu2,
|
||||
__clc_fma(-u2, vu1, __clc_fma(-u1, vu1, v))),
|
||||
u));
|
||||
}
|
||||
|
||||
q3 = vbyu < 0x1.d12ed0af1a27fp-27 ? q3 : q5;
|
||||
q4 = vbyu < 0x1.d12ed0af1a27fp-27 ? q4 : q6;
|
||||
|
||||
q1 = vbyu > 0.0625 ? q1 : q3;
|
||||
q2 = vbyu > 0.0625 ? q2 : q4;
|
||||
|
||||
// Tidy-up according to which quadrant the arguments lie in
|
||||
__CLC_GENTYPE res1, res2, res3, res4;
|
||||
q1 = swap_vu ? piby2_head - q1 : q1;
|
||||
q2 = swap_vu ? piby2_tail - q2 : q2;
|
||||
q1 = xneg ? pi_head - q1 : q1;
|
||||
q2 = xneg ? pi_tail - q2 : q2;
|
||||
q1 = MATH_DIVIDE(q1 + q2, pi);
|
||||
res4 = yneg ? -q1 : q1;
|
||||
|
||||
res1 = yneg ? -0.75 : 0.75;
|
||||
res2 = yneg ? -0.25 : 0.25;
|
||||
res3 = xneg ? res1 : res2;
|
||||
|
||||
res3 = __clc_select(res4, res3,
|
||||
__CLC_CONVERT_LONGN(__clc_isinf(y2) & __clc_isinf(x2)));
|
||||
res1 = yneg ? -1.0 : 1.0;
|
||||
|
||||
// abs(x)/abs(y) > 2^56 and x < 0
|
||||
res3 = diffexp < -56 && xneg ? res1 : res3;
|
||||
|
||||
res4 = MATH_DIVIDE(MATH_DIVIDE(y, x), pi);
|
||||
// x positive and dominant over y by a factor of 2^28
|
||||
res3 = diffexp < -28 && xneg == 0 ? res4 : res3;
|
||||
|
||||
// abs(y)/abs(x) > 2^56
|
||||
res4 = yneg ? -0.5 : 0.5; // atan(y/x) is insignificant compared to piby2
|
||||
res3 = diffexp > 56 ? res4 : res3;
|
||||
|
||||
res3 = x2 == 0.0 ? res4 : res3; // Zero x gives +- pi/2 depending on sign of y
|
||||
res4 = xneg ? res1 : y2;
|
||||
|
||||
// Zero y gives +-0 for positive x and +-pi for negative x
|
||||
res3 = y2 == 0.0 ? res4 : res3;
|
||||
res3 = __clc_isnan(y2) ? y2 : res3;
|
||||
res3 = __clc_isnan(x2) ? x2 : res3;
|
||||
|
||||
return res3;
|
||||
return __clc_copysign(a, y);
|
||||
}
|
||||
|
||||
#elif __CLC_FPSIZE == 16
|
||||
|
||||
_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_atan2pi(__CLC_GENTYPE x,
|
||||
__CLC_GENTYPE y) {
|
||||
return __CLC_CONVERT_GENTYPE(
|
||||
__clc_atan2pi(__CLC_CONVERT_FLOATN(x), __CLC_CONVERT_FLOATN(y)));
|
||||
_CLC_DEF _CLC_OVERLOAD _CLC_CONST __CLC_HALFN __clc_atan2pi(__CLC_HALFN y,
|
||||
__CLC_HALFN x) {
|
||||
__CLC_HALFN ax = __clc_fabs(x);
|
||||
__CLC_HALFN ay = __clc_fabs(y);
|
||||
__CLC_HALFN v = __clc_fmin(ax, ay);
|
||||
__CLC_HALFN u = __clc_fmax(ax, ay);
|
||||
|
||||
__CLC_HALFN vbyu = v / u;
|
||||
|
||||
__CLC_HALFN a = __clc_atanpi_reduced(vbyu);
|
||||
|
||||
__CLC_HALFN at = 0.5h - a;
|
||||
a = ay > ax ? at : a;
|
||||
at = 1.0h - a;
|
||||
a = x < 0.0h ? at : a;
|
||||
|
||||
at = __clc_signbit(x) ? 1.0h : 0.0h;
|
||||
a = y == 0.0h ? at : a;
|
||||
|
||||
// x and y are +- Inf
|
||||
at = x < 0.0h ? 0.75h : 0.25h;
|
||||
a = (__clc_isinf(x) & __clc_isinf(y)) ? at : a;
|
||||
|
||||
// x or y is NaN
|
||||
a = __clc_isunordered(x, y) ? HALF_NAN : a;
|
||||
|
||||
return __clc_copysign(a, y);
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user