[libclc] Move exp10 to the CLC library

libclc/clc/include/clc/math/clc_exp10.h

@@ -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_EXP10_H__
+#define __CLC_MATH_CLC_EXP10_H__
+
+#define __CLC_BODY <clc/math/unary_decl.inc>
+#define __CLC_FUNCTION __clc_exp10
+
+#include <clc/math/gentype.inc>
+
+#undef __CLC_BODY
+#undef __CLC_FUNCTION
+
+#endif // __CLC_MATH_CLC_EXP10_H__

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

@@ -64,13 +64,13 @@ TABLE_FUNCTION_DECL(float2, log10_tbl);
 TABLE_FUNCTION_DECL(uint4, pibits_tbl);
 TABLE_FUNCTION_DECL(float2, sinhcosh_tbl);
 TABLE_FUNCTION_DECL(float2, cbrt_tbl);
-TABLE_FUNCTION_DECL(float, exp_tbl);
 
 CLC_TABLE_FUNCTION_DECL(float, log_inv_tbl_ep_head);
 CLC_TABLE_FUNCTION_DECL(float, log_inv_tbl_ep_tail);
 CLC_TABLE_FUNCTION_DECL(float, loge_tbl_lo);
 CLC_TABLE_FUNCTION_DECL(float, loge_tbl_hi);
 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);
 

libclc/clc/lib/generic/SOURCES

@@ -32,6 +32,7 @@ math/clc_ceil.cl
 math/clc_copysign.cl
 math/clc_cospi.cl
 math/clc_ep_log.cl
+math/clc_exp10.cl
 math/clc_fabs.cl
 math/clc_fma.cl
 math/clc_fmod.cl

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

@@ -0,0 +1,21 @@
+//===----------------------------------------------------------------------===//
+//
+// 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_fma.h>
+#include <clc/math/clc_ldexp.h>
+#include <clc/math/clc_mad.h>
+#include <clc/math/clc_subnormal_config.h>
+#include <clc/math/math.h>
+#include <clc/math/tables.h>
+#include <clc/relational/clc_isnan.h>
+
+#define __CLC_BODY <clc_exp10.inc>
+#include <clc/math/gentype.inc>

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

@@ -0,0 +1,155 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+//
+//  Algorithm:
+//
+//  e^x = 2^(x/ln(2)) = 2^(x*(64/ln(2))/64)
+//
+//  x*(64/ln(2)) = n + f, |f| <= 0.5, n is integer
+//  n = 64*m + j,   0 <= j < 64
+//
+//  e^x = 2^((64*m + j + f)/64)
+//      = (2^m) * (2^(j/64)) * 2^(f/64)
+//      = (2^m) * (2^(j/64)) * e^(f*(ln(2)/64))
+//
+//  f = x*(64/ln(2)) - n
+//  r = f*(ln(2)/64) = x - n*(ln(2)/64)
+//
+//  e^x = (2^m) * (2^(j/64)) * e^r
+//
+//  (2^(j/64)) is precomputed
+//
+//  e^r = 1 + r + (r^2)/2! + (r^3)/3! + (r^4)/4! + (r^5)/5!
+//  e^r = 1 + q
+//
+//  q = r + (r^2)/2! + (r^3)/3! + (r^4)/4! + (r^5)/5!
+//
+//  e^x = (2^m) * ( (2^(j/64)) + q*(2^(j/64)) )
+//
+//===----------------------------------------------------------------------===//
+
+#if __CLC_FPSIZE == 32
+
+_CLC_DEF _CLC_OVERLOAD __CLC_GENTYPE __clc_exp10(__CLC_GENTYPE x) {
+  // 128*log2/log10 : 38.53183944498959
+  const __CLC_GENTYPE X_MAX = 0x1.344134p+5f;
+  // -149*log2/log10 : -44.8534693539332
+  const __CLC_GENTYPE X_MIN = -0x1.66d3e8p+5f;
+  // 64*log10/log2 : 212.6033980727912
+  const __CLC_GENTYPE R_64_BY_LOG10_2 = 0x1.a934f0p+7f;
+  // log2/(64 * log10) lead : 0.004699707
+  const __CLC_GENTYPE R_LOG10_2_BY_64_LD = 0x1.340000p-8f;
+  // log2/(64 * log10) tail : 0.00000388665057
+  const __CLC_GENTYPE R_LOG10_2_BY_64_TL = 0x1.04d426p-18f;
+  const __CLC_GENTYPE R_LN10 = 0x1.26bb1cp+1f;
+
+  __CLC_INTN return_nan = __clc_isnan(x);
+  __CLC_INTN return_inf = x > X_MAX;
+  __CLC_INTN return_zero = x < X_MIN;
+
+  __CLC_INTN n = __CLC_CONVERT_INTN(x * R_64_BY_LOG10_2);
+
+  __CLC_GENTYPE fn = __CLC_CONVERT_GENTYPE(n);
+  __CLC_INTN j = n & 0x3f;
+  __CLC_INTN m = n >> 6;
+  __CLC_INTN m2 = m << EXPSHIFTBITS_SP32;
+  __CLC_GENTYPE r;
+
+  r = R_LN10 *
+      __clc_mad(fn, -R_LOG10_2_BY_64_TL, __clc_mad(fn, -R_LOG10_2_BY_64_LD, x));
+
+  // Truncated Taylor series for e^r
+  __CLC_GENTYPE z2 =
+      __clc_mad(__clc_mad(__clc_mad(r, 0x1.555556p-5f, 0x1.555556p-3f), r,
+                          0x1.000000p-1f),
+                r * r, r);
+
+  __CLC_GENTYPE two_to_jby64 = USE_TABLE(exp_tbl, j);
+  z2 = __clc_mad(two_to_jby64, z2, two_to_jby64);
+
+  __CLC_GENTYPE z2s = z2 * __CLC_AS_GENTYPE((__CLC_UINTN)0x1 << (m + 149));
+  __CLC_GENTYPE z2n = __CLC_AS_GENTYPE(__CLC_AS_INTN(z2) + m2);
+  z2 = m <= -126 ? z2s : z2n;
+
+  z2 = return_inf ? __CLC_AS_GENTYPE((__CLC_UINTN)PINFBITPATT_SP32) : z2;
+  z2 = return_zero ? 0.0f : z2;
+  z2 = return_nan ? x : z2;
+  return z2;
+}
+
+#elif __CLC_FPSIZE == 64
+
+_CLC_DEF _CLC_OVERLOAD __CLC_GENTYPE __clc_exp10(__CLC_GENTYPE x) {
+  // 1024*ln(2)/ln(10)
+  const __CLC_GENTYPE X_MAX = 0x1.34413509f79ffp+8;
+  // -1074*ln(2)/ln(10)
+  const __CLC_GENTYPE X_MIN = -0x1.434e6420f4374p+8;
+  // 64*ln(10)/ln(2)
+  const __CLC_GENTYPE R_64_BY_LOG10_2 = 0x1.a934f0979a371p+7;
+  // head ln(2)/(64*ln(10))
+  const __CLC_GENTYPE R_LOG10_2_BY_64_LD = 0x1.3441350000000p-8;
+  // tail ln(2)/(64*ln(10))
+  const __CLC_GENTYPE R_LOG10_2_BY_64_TL = 0x1.3ef3fde623e25p-37;
+  // ln(10)
+  const __CLC_GENTYPE R_LN10 = 0x1.26bb1bbb55516p+1;
+
+  __CLC_INTN n = __CLC_CONVERT_INTN(x * R_64_BY_LOG10_2);
+
+  __CLC_GENTYPE dn = __CLC_CONVERT_GENTYPE(n);
+
+  __CLC_INTN j = n & 0x3f;
+  __CLC_INTN m = n >> 6;
+
+  __CLC_GENTYPE r = R_LN10 * __clc_fma(-R_LOG10_2_BY_64_TL, dn,
+                                       __clc_fma(-R_LOG10_2_BY_64_LD, dn, x));
+
+  // 6 term tail of Taylor expansion of e^r
+  __CLC_GENTYPE z2 =
+      r * __clc_fma(
+              r,
+              __clc_fma(r,
+                        __clc_fma(r,
+                                  __clc_fma(r,
+                                            __clc_fma(r, 0x1.6c16c16c16c17p-10,
+                                                      0x1.1111111111111p-7),
+                                            0x1.5555555555555p-5),
+                                  0x1.5555555555555p-3),
+                        0x1.0000000000000p-1),
+              1.0);
+
+  __CLC_GENTYPE tv0 = USE_TABLE(two_to_jby64_ep_tbl_head, j);
+  __CLC_GENTYPE tv1 = USE_TABLE(two_to_jby64_ep_tbl_tail, j);
+  z2 = __clc_fma(tv0 + tv1, z2, tv1) + tv0;
+
+  __CLC_INTN small_value =
+      (m < -1022) || ((m == -1022) && __CLC_CONVERT_INTN(z2 < 1.0));
+
+  __CLC_INTN n1 = m >> 2;
+  __CLC_INTN n2 = m - n1;
+  __CLC_GENTYPE z3 =
+      z2 * __CLC_AS_GENTYPE((__CLC_CONVERT_LONGN(n1) + 1023) << 52);
+  z3 *= __CLC_AS_GENTYPE((__CLC_CONVERT_LONGN(n2) + 1023) << 52);
+
+  z2 = __clc_ldexp(z2, m);
+  z2 = __CLC_CONVERT_LONGN(small_value) ? z3 : z2;
+
+  z2 = __clc_isnan(x) ? x : z2;
+
+  z2 = x > X_MAX ? __CLC_AS_GENTYPE((__CLC_ULONGN)PINFBITPATT_DP64) : z2;
+  z2 = x < X_MIN ? 0.0 : z2;
+
+  return z2;
+}
+
+#elif __CLC_FPSIZE == 16
+
+_CLC_OVERLOAD _CLC_DEF __CLC_GENTYPE __clc_exp10(__CLC_GENTYPE x) {
+  return __CLC_CONVERT_GENTYPE(__clc_exp10(__CLC_CONVERT_FLOATN(x)));
+}
+
+#endif

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

@@ -197,6 +197,28 @@ DECLARE_TABLE(float, LOG_INV_TBL, 129) = {
 
 CLC_TABLE_FUNCTION(float, LOG_INV_TBL, log_inv_tbl);
 
+DECLARE_TABLE(float, EXP_TBL, 65) = {
+    0x1.000000p+0f, 0x1.02c9a4p+0f, 0x1.059b0ep+0f, 0x1.087452p+0f,
+    0x1.0b5586p+0f, 0x1.0e3ec4p+0f, 0x1.11301ep+0f, 0x1.1429aap+0f,
+    0x1.172b84p+0f, 0x1.1a35bep+0f, 0x1.1d4874p+0f, 0x1.2063b8p+0f,
+    0x1.2387a6p+0f, 0x1.26b456p+0f, 0x1.29e9e0p+0f, 0x1.2d285ap+0f,
+    0x1.306fe0p+0f, 0x1.33c08cp+0f, 0x1.371a74p+0f, 0x1.3a7db4p+0f,
+    0x1.3dea64p+0f, 0x1.4160a2p+0f, 0x1.44e086p+0f, 0x1.486a2cp+0f,
+    0x1.4bfdaep+0f, 0x1.4f9b28p+0f, 0x1.5342b6p+0f, 0x1.56f474p+0f,
+    0x1.5ab07ep+0f, 0x1.5e76f2p+0f, 0x1.6247ecp+0f, 0x1.662388p+0f,
+    0x1.6a09e6p+0f, 0x1.6dfb24p+0f, 0x1.71f75ep+0f, 0x1.75feb6p+0f,
+    0x1.7a1148p+0f, 0x1.7e2f34p+0f, 0x1.82589ap+0f, 0x1.868d9ap+0f,
+    0x1.8ace54p+0f, 0x1.8f1aeap+0f, 0x1.93737cp+0f, 0x1.97d82ap+0f,
+    0x1.9c4918p+0f, 0x1.a0c668p+0f, 0x1.a5503cp+0f, 0x1.a9e6b6p+0f,
+    0x1.ae89fap+0f, 0x1.b33a2cp+0f, 0x1.b7f770p+0f, 0x1.bcc1eap+0f,
+    0x1.c199bep+0f, 0x1.c67f12p+0f, 0x1.cb720ep+0f, 0x1.d072d4p+0f,
+    0x1.d5818ep+0f, 0x1.da9e60p+0f, 0x1.dfc974p+0f, 0x1.e502eep+0f,
+    0x1.ea4afap+0f, 0x1.efa1bep+0f, 0x1.f50766p+0f, 0x1.fa7c18p+0f,
+    0x1.000000p+1f,
+};
+
+CLC_TABLE_FUNCTION(float, EXP_TBL, exp_tbl);
+
 DECLARE_TABLE(float, EXP_TBL_EP_HEAD, 65) = {
     0x1.000000p+0f, 0x1.02c000p+0f, 0x1.058000p+0f, 0x1.084000p+0f,
     0x1.0b4000p+0f, 0x1.0e0000p+0f, 0x1.110000p+0f, 0x1.140000p+0f,

libclc/clspv/lib/SOURCES

@@ -16,7 +16,6 @@ subnormal_config.cl
 ../../generic/lib/math/atanh.cl
 ../../generic/lib/math/atanpi.cl
 ../../generic/lib/math/cbrt.cl
-../../generic/lib/math/clc_exp10.cl
 ../../generic/lib/math/clc_tan.cl
 ../../generic/lib/math/cos.cl
 ../../generic/lib/math/cosh.cl

libclc/generic/lib/SOURCES

@@ -100,7 +100,6 @@ math/exp.cl
 math/exp_helper.cl
 math/expm1.cl
 math/exp2.cl
-math/clc_exp10.cl
 math/exp10.cl
 math/fabs.cl
 math/fdim.cl