[libclc] Move cbrt to the CLC library; vectorize

libclc/clc/include/clc/math/clc_cbrt.inc

@@ -0,0 +1,20 @@
+//===----------------------------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef __CLC_MATH_CLC_CBRT_H__
+#define __CLC_MATH_CLC_CBRT_H__
+
+#define __CLC_BODY <clc/math/unary_decl.inc>
+#define __CLC_FUNCTION __clc_cbrt
+
+#include <clc/math/gentype.inc>
+
+#undef __CLC_BODY
+#undef __CLC_FUNCTION
+
+#endif // __CLC_MATH_CLC_CBRT_H__

libclc/clc/include/clc/math/tables.h

@@ -63,7 +63,6 @@ TABLE_FUNCTION_DECL(float2, log2_tbl);
 TABLE_FUNCTION_DECL(float2, log10_tbl);
 TABLE_FUNCTION_DECL(uint4, pibits_tbl);
 TABLE_FUNCTION_DECL(float2, sinhcosh_tbl);
-TABLE_FUNCTION_DECL(float2, cbrt_tbl);
 
 CLC_TABLE_FUNCTION_DECL(float, log_inv_tbl_ep_head);
 CLC_TABLE_FUNCTION_DECL(float, log_inv_tbl_ep_tail);
@@ -73,6 +72,8 @@ CLC_TABLE_FUNCTION_DECL(float, log_inv_tbl);
 CLC_TABLE_FUNCTION_DECL(float, exp_tbl);
 CLC_TABLE_FUNCTION_DECL(float, exp_tbl_ep_head);
 CLC_TABLE_FUNCTION_DECL(float, exp_tbl_ep_tail);
+CLC_TABLE_FUNCTION_DECL(float, cbrt_tbl_head);
+CLC_TABLE_FUNCTION_DECL(float, cbrt_tbl_tail);
 
 #ifdef cl_khr_fp64
 
@@ -86,9 +87,11 @@ CLC_TABLE_FUNCTION_DECL(double, two_to_jby64_ep_tbl_head);
 CLC_TABLE_FUNCTION_DECL(double, two_to_jby64_ep_tbl_tail);
 TABLE_FUNCTION_DECL(double2, sinh_tbl);
 TABLE_FUNCTION_DECL(double2, cosh_tbl);
-TABLE_FUNCTION_DECL(double, cbrt_inv_tbl);
-TABLE_FUNCTION_DECL(double2, cbrt_dbl_tbl);
-TABLE_FUNCTION_DECL(double2, cbrt_rem_tbl);
+CLC_TABLE_FUNCTION_DECL(double, cbrt_inv_tbl);
+CLC_TABLE_FUNCTION_DECL(double, cbrt_dbl_tbl_head);
+CLC_TABLE_FUNCTION_DECL(double, cbrt_dbl_tbl_tail);
+CLC_TABLE_FUNCTION_DECL(double, cbrt_rem_tbl_head);
+CLC_TABLE_FUNCTION_DECL(double, cbrt_rem_tbl_tail);
 CLC_TABLE_FUNCTION_DECL(double, powlog_tbl_head);
 CLC_TABLE_FUNCTION_DECL(double, powlog_tbl_tail);
 CLC_TABLE_FUNCTION_DECL(double, log_f_inv_tbl_head);

libclc/clc/lib/generic/SOURCES

@@ -28,6 +28,7 @@ math/clc_atan2.cl
 math/clc_atan2pi.cl
 math/clc_atanh.cl
 math/clc_atanpi.cl
+math/clc_cbrt.cl
 math/clc_ceil.cl
 math/clc_copysign.cl
 math/clc_cospi.cl

libclc/clc/lib/generic/math/clc_cbrt.cl

@@ -0,0 +1,22 @@
+//===----------------------------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include <clc/clc_convert.h>
+#include <clc/clcmacro.h>
+#include <clc/internal/clc.h>
+#include <clc/math/clc_copysign.h>
+#include <clc/math/clc_fabs.h>
+#include <clc/math/clc_fma.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>
+
+#define __CLC_BODY <clc_cbrt.inc>
+#include <clc/math/gentype.inc>

libclc/clc/lib/generic/math/clc_cbrt.inc

@@ -0,0 +1,135 @@
+//===----------------------------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#if __CLC_FPSIZE == 32
+
+_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_cbrt(__CLC_GENTYPE x) {
+  __CLC_UINTN xi = __CLC_AS_UINTN(x);
+  __CLC_UINTN axi = __CLC_AS_UINTN(__clc_fabs(x));
+  __CLC_UINTN xsign = axi ^ xi;
+  xi = axi;
+
+  __CLC_INTN m = __CLC_AS_INTN(xi >> EXPSHIFTBITS_SP32) - EXPBIAS_SP32;
+
+  // Treat subnormals
+  __CLC_UINTN xisub = __CLC_AS_UINTN(__CLC_AS_GENTYPE(xi | 0x3f800000) - 1.0f);
+  __CLC_INTN msub = __CLC_AS_INTN(xisub >> EXPSHIFTBITS_SP32) - 253;
+  __CLC_INTN c = m == -127;
+  xi = c ? xisub : xi;
+  m = c ? msub : m;
+
+  __CLC_INTN m3 = m / 3;
+  __CLC_INTN rem = m - m3 * 3;
+  __CLC_GENTYPE mf = __CLC_AS_GENTYPE((m3 + EXPBIAS_SP32) << EXPSHIFTBITS_SP32);
+
+  __CLC_UINTN indx = (xi & 0x007f0000) + ((xi & 0x00008000) << 1);
+  __CLC_GENTYPE f = __CLC_AS_GENTYPE((xi & MANTBITS_SP32) | 0x3f000000) -
+                    __CLC_AS_GENTYPE(indx | 0x3f000000);
+
+  indx >>= 16;
+  __CLC_GENTYPE r = f * USE_TABLE(log_inv_tbl, __CLC_AS_INTN(indx));
+  __CLC_GENTYPE poly = __clc_mad(__clc_mad(r, 0x1.f9add4p-5f, -0x1.c71c72p-4f),
+                                 r * r, r * 0x1.555556p-2f);
+
+  // This could also be done with a 5-element table
+  __CLC_GENTYPE remH = 0x1.428000p-1f;
+  __CLC_GENTYPE remT = 0x1.45f31ap-14f;
+
+  remH = rem == -1 ? 0x1.964000p-1f : remH;
+  remT = rem == -1 ? 0x1.fea53ep-13f : remT;
+
+  remH = rem == 0 ? 0x1.000000p+0f : remH;
+  remT = rem == 0 ? 0x0.000000p+0f : remT;
+
+  remH = rem == 1 ? 0x1.428000p+0f : remH;
+  remT = rem == 1 ? 0x1.45f31ap-13f : remT;
+
+  remH = rem == 2 ? 0x1.964000p+0f : remH;
+  remT = rem == 2 ? 0x1.fea53ep-12f : remT;
+
+  __CLC_GENTYPE cbrtH = USE_TABLE(cbrt_tbl_head, __CLC_AS_INTN(indx));
+  __CLC_GENTYPE cbrtT = USE_TABLE(cbrt_tbl_tail, __CLC_AS_INTN(indx));
+
+  __CLC_GENTYPE bH = cbrtH * remH;
+  __CLC_GENTYPE bT =
+      __clc_mad(cbrtH, remT, __clc_mad(cbrtT, remH, cbrtT * remT));
+
+  __CLC_GENTYPE z = __clc_mad(poly, bH, __clc_mad(poly, bT, bT)) + bH;
+  z *= mf;
+  z = __CLC_AS_GENTYPE(__CLC_AS_UINTN(z) | xsign);
+  c = axi >= EXPBITS_SP32 || axi == 0;
+  z = c ? x : z;
+  return z;
+}
+
+#elif __CLC_FPSIZE == 64
+
+_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_cbrt(__CLC_GENTYPE x) {
+  __CLC_LONGN return_x = __clc_isinf(x) || __clc_isnan(x) || x == 0.0;
+  __CLC_ULONGN ux = __CLC_AS_ULONGN(__clc_fabs(x));
+  __CLC_INTN m = __CLC_CONVERT_INTN(ux >> EXPSHIFTBITS_DP64) - 1023;
+
+  // Treat subnormals
+  __CLC_ULONGN uxs =
+      __CLC_AS_ULONGN(__CLC_AS_GENTYPE(0x3ff0000000000000UL | ux) - 1.0);
+  __CLC_INTN ms = m + __CLC_CONVERT_INTN(uxs >> EXPSHIFTBITS_DP64) - 1022;
+
+  __CLC_INTN c = m == -1023;
+  ux = __CLC_CONVERT_LONGN(c) ? uxs : ux;
+  m = c ? ms : m;
+
+  __CLC_INTN mby3 = m / 3;
+  __CLC_INTN rem = m - 3 * mby3;
+
+  __CLC_GENTYPE mf = __CLC_AS_GENTYPE(__CLC_CONVERT_ULONGN(mby3 + 1023) << 52);
+
+  ux &= 0x000fffffffffffffUL;
+  __CLC_GENTYPE Y = __CLC_AS_GENTYPE(0x3fe0000000000000UL | ux);
+
+  // nearest integer
+  __CLC_INTN index = __CLC_CONVERT_INTN(ux >> 43);
+  index = (0x100 | (index >> 1)) + (index & 1);
+  __CLC_GENTYPE F = __CLC_CONVERT_GENTYPE(index) * 0x1.0p-9;
+
+  __CLC_GENTYPE f = Y - F;
+  __CLC_GENTYPE r = f * USE_TABLE(cbrt_inv_tbl, index - 256);
+
+  __CLC_GENTYPE z =
+      r * __clc_fma(
+              r,
+              __clc_fma(r,
+                        __clc_fma(r,
+                                  __clc_fma(r,
+                                            __clc_fma(r, -0x1.8090d6221a247p-6,
+                                                      0x1.ee7113506ac13p-6),
+                                            -0x1.511e8d2b3183bp-5),
+                                  0x1.f9add3c0ca458p-5),
+                        -0x1.c71c71c71c71cp-4),
+              0x1.5555555555555p-2);
+
+  __CLC_GENTYPE Rem_h = USE_TABLE(cbrt_rem_tbl_head, rem + 2);
+  __CLC_GENTYPE Rem_t = USE_TABLE(cbrt_rem_tbl_tail, rem + 2);
+
+  __CLC_GENTYPE F_h = USE_TABLE(cbrt_dbl_tbl_head, index - 256);
+  __CLC_GENTYPE F_t = USE_TABLE(cbrt_dbl_tbl_tail, index - 256);
+
+  __CLC_GENTYPE b_h = F_h * Rem_h;
+  __CLC_GENTYPE b_t = __clc_fma(Rem_t, F_h, __clc_fma(F_t, Rem_h, F_t * Rem_t));
+
+  __CLC_GENTYPE ans = __clc_fma(z, b_h, __clc_fma(z, b_t, b_t)) + b_h;
+  ans = __clc_copysign(ans * mf, x);
+  return return_x ? x : ans;
+}
+
+#elif __CLC_FPSIZE == 16
+
+_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_cbrt(__CLC_GENTYPE x) {
+  return __CLC_CONVERT_GENTYPE(__clc_cbrt(__CLC_CONVERT_FLOATN(x)));
+}
+
+#endif