[Libclc-dev] [PATCH 1/2] math: Add expm1 builtin function
Aaron Watry via Libclc-dev
libclc-dev at lists.llvm.org
Mon Jan 16 18:13:22 PST 2017
Ported from the amd-builtins branch.
Signed-off-by: Aaron Watry <awatry at gmail.com>
CC: Tom Stellard <thomas.stellard at amd.com>
---
generic/include/clc/clc.h | 1 +
generic/include/clc/math/expm1.h | 9 +++
generic/lib/SOURCES | 1 +
generic/lib/math/expm1.cl | 142 +++++++++++++++++++++++++++++++++++++++
generic/lib/math/tables.cl | 138 +++++++++++++++++++++++++++++++++++++
generic/lib/math/tables.h | 2 +
6 files changed, 293 insertions(+)
create mode 100644 generic/include/clc/math/expm1.h
create mode 100644 generic/lib/math/expm1.cl
diff --git a/generic/include/clc/clc.h b/generic/include/clc/clc.h
index 5bdb09f..f95d7a7 100644
--- a/generic/include/clc/clc.h
+++ b/generic/include/clc/clc.h
@@ -53,6 +53,7 @@
#include <clc/math/erf.h>
#include <clc/math/erfc.h>
#include <clc/math/exp.h>
+#include <clc/math/expm1.h>
#include <clc/math/exp10.h>
#include <clc/math/exp2.h>
#include <clc/math/fabs.h>
diff --git a/generic/include/clc/math/expm1.h b/generic/include/clc/math/expm1.h
new file mode 100644
index 0000000..522ddfe
--- /dev/null
+++ b/generic/include/clc/math/expm1.h
@@ -0,0 +1,9 @@
+#undef exp
+
+#define __CLC_BODY <clc/math/unary_decl.inc>
+#define __CLC_FUNCTION expm1
+
+#include <clc/math/gentype.inc>
+
+#undef __CLC_BODY
+#undef __CLC_FUNCTION
diff --git a/generic/lib/SOURCES b/generic/lib/SOURCES
index ecd2e73..74bea25 100644
--- a/generic/lib/SOURCES
+++ b/generic/lib/SOURCES
@@ -83,6 +83,7 @@ math/erf.cl
math/erfc.cl
math/exp.cl
math/exp_helper.cl
+math/expm1.cl
math/exp2.cl
math/exp10.cl
math/fdim.cl
diff --git a/generic/lib/math/expm1.cl b/generic/lib/math/expm1.cl
new file mode 100644
index 0000000..9a3a907
--- /dev/null
+++ b/generic/lib/math/expm1.cl
@@ -0,0 +1,142 @@
+#include <clc/clc.h>
+
+#include "math.h"
+#include "tables.h"
+#include "../clcmacro.h"
+
+/* Refer to the exp routine for the underlying algorithm */
+
+_CLC_OVERLOAD _CLC_DEF float expm1(float x) {
+ const float X_MAX = 0x1.62e42ep+6f; // 128*log2 : 88.722839111673
+ const float X_MIN = -0x1.9d1da0p+6f; // -149*log2 : -103.27892990343184
+
+ const float R_64_BY_LOG2 = 0x1.715476p+6f; // 64/log2 : 92.332482616893657
+ const float R_LOG2_BY_64_LD = 0x1.620000p-7f; // log2/64 lead: 0.0108032227
+ const float R_LOG2_BY_64_TL = 0x1.c85fdep-16f; // log2/64 tail: 0.0000272020388
+
+ uint xi = as_uint(x);
+ int n = (int)(x * R_64_BY_LOG2);
+ float fn = (float)n;
+
+ int j = n & 0x3f;
+ int m = n >> 6;
+
+ float r = mad(fn, -R_LOG2_BY_64_TL, mad(fn, -R_LOG2_BY_64_LD, x));
+
+ // Truncated Taylor series
+ float z2 = mad(r*r, mad(r, mad(r, 0x1.555556p-5f, 0x1.555556p-3f), 0.5f), r);
+
+ float m2 = as_float((m + EXPBIAS_SP32) << EXPSHIFTBITS_SP32);
+ float2 tv = USE_TABLE(exp_tbl_ep, j);
+
+ float two_to_jby64_h = tv.s0 * m2;
+ float two_to_jby64_t = tv.s1 * m2;
+ float two_to_jby64 = two_to_jby64_h + two_to_jby64_t;
+
+ z2 = mad(z2, two_to_jby64, two_to_jby64_t) + (two_to_jby64_h - 1.0f);
+ //Make subnormals work
+ z2 = x == 0.f ? x : z2;
+ z2 = x < X_MIN | m < -24 ? -1.0f : z2;
+ z2 = x > X_MAX ? as_float(PINFBITPATT_SP32) : z2;
+ z2 = isnan(x) ? x : z2;
+
+ return z2;
+}
+
+_CLC_UNARY_VECTORIZE(_CLC_OVERLOAD _CLC_DEF, float, expm1, float)
+
+#ifdef cl_khr_fp64
+
+#include "exp_helper.h"
+
+#pragma OPENCL EXTENSION cl_khr_fp64 : enable
+
+_CLC_OVERLOAD _CLC_DEF double expm1(double x) {
+ const double max_expm1_arg = 709.8;
+ const double min_expm1_arg = -37.42994775023704;
+ const double log_OnePlus_OneByFour = 0.22314355131420976; //0x3FCC8FF7C79A9A22 = log(1+1/4)
+ const double log_OneMinus_OneByFour = -0.28768207245178096; //0xBFD269621134DB93 = log(1-1/4)
+ const double sixtyfour_by_lnof2 = 92.33248261689366; //0x40571547652b82fe
+ const double lnof2_by_64_head = 0.010830424696223417; //0x3f862e42fefa0000
+ const double lnof2_by_64_tail = 2.5728046223276688e-14; //0x3d1cf79abc9e3b39
+
+ // First, assume log(1-1/4) < x < log(1+1/4) i.e -0.28768 < x < 0.22314
+ double u = as_double(as_ulong(x) & 0xffffffffff000000UL);
+ double v = x - u;
+ double y = u * u * 0.5;
+ double z = v * (x + u) * 0.5;
+
+ double q = fma(x,
+ fma(x,
+ fma(x,
+ fma(x,
+ fma(x,
+ fma(x,
+ fma(x,
+ fma(x,2.4360682937111612e-8, 2.7582184028154370e-7),
+ 2.7558212415361945e-6),
+ 2.4801576918453420e-5),
+ 1.9841269447671544e-4),
+ 1.3888888890687830e-3),
+ 8.3333333334012270e-3),
+ 4.1666666666665560e-2),
+ 1.6666666666666632e-1);
+ q *= x * x * x;
+
+ double z1g = (u + y) + (q + (v + z));
+ double z1 = x + (y + (q + z));
+ z1 = y >= 0x1.0p-7 ? z1g : z1;
+
+ // Now assume outside interval around 0
+ int n = (int)(x * sixtyfour_by_lnof2);
+ int j = n & 0x3f;
+ int m = n >> 6;
+
+ double2 tv = USE_TABLE(two_to_jby64_ep_tbl, j);
+ double f1 = tv.s0;
+ double f2 = tv.s1;
+ double f = f1 + f2;
+
+ double dn = -n;
+ double r = fma(dn, lnof2_by_64_tail, fma(dn, lnof2_by_64_head, x));
+
+ q = fma(r,
+ fma(r,
+ fma(r,
+ fma(r, 1.38889490863777199667e-03, 8.33336798434219616221e-03),
+ 4.16666666662260795726e-02),
+ 1.66666666665260878863e-01),
+ 5.00000000000000008883e-01);
+ q = fma(r*r, q, r);
+
+ double twopm = as_double((long)(m + EXPBIAS_DP64) << EXPSHIFTBITS_DP64);
+ double twopmm = as_double((long)(EXPBIAS_DP64 - m) << EXPSHIFTBITS_DP64);
+
+ // Computations for m > 52, including where result is close to Inf
+ ulong uval = as_ulong(0x1.0p+1023 * (f1 + (f * q + (f2))));
+ int e = (int)(uval >> EXPSHIFTBITS_DP64) + 1;
+
+ double zme1024 = as_double(((long)e << EXPSHIFTBITS_DP64) | (uval & MANTBITS_DP64));
+ zme1024 = e == 2047 ? as_double(PINFBITPATT_DP64) : zme1024;
+
+ double zmg52 = twopm * (f1 + fma(f, q, f2 - twopmm));
+ zmg52 = m == 1024 ? zme1024 : zmg52;
+
+ // For m < 53
+ double zml53 = twopm * ((f1 - twopmm) + fma(f1, q, f2*(1.0 + q)));
+
+ // For m < -7
+ double zmln7 = fma(twopm, f1 + fma(f, q, f2), -1.0);
+
+ z = m < 53 ? zml53 : zmg52;
+ z = m < -7 ? zmln7 : z;
+ z = x > log_OneMinus_OneByFour & x < log_OnePlus_OneByFour ? z1 : z;
+ z = x > max_expm1_arg ? as_double(PINFBITPATT_DP64) : z;
+ z = x < min_expm1_arg ? -1.0 : z;
+
+ return z;
+}
+
+_CLC_UNARY_VECTORIZE(_CLC_OVERLOAD _CLC_DEF, double, expm1, double)
+
+#endif
diff --git a/generic/lib/math/tables.cl b/generic/lib/math/tables.cl
index c1dd85c..3e1e4db 100644
--- a/generic/lib/math/tables.cl
+++ b/generic/lib/math/tables.cl
@@ -608,6 +608,142 @@ DECLARE_TABLE(float2, CBRT_TBL, 129) = {
(float2)(0x1.428000p+0f, 0x1.45f31ap-13f)
};
+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,
+};
+
+DECLARE_TABLE(float2, EXP_TBL_EP, 65) = {
+ (float2) (0x1.000000p+0f, 0x0.000000p+0f),
+ (float2) (0x1.02c000p+0f, 0x1.347ceep-13f),
+ (float2) (0x1.058000p+0f, 0x1.b0d314p-12f),
+ (float2) (0x1.084000p+0f, 0x1.a28c3ap-11f),
+ (float2) (0x1.0b4000p+0f, 0x1.586cf8p-12f),
+ (float2) (0x1.0e0000p+0f, 0x1.f61968p-11f),
+ (float2) (0x1.110000p+0f, 0x1.80e808p-11f),
+ (float2) (0x1.140000p+0f, 0x1.4d5754p-11f),
+ (float2) (0x1.170000p+0f, 0x1.5c1e3ep-11f),
+ (float2) (0x1.1a0000p+0f, 0x1.adf5b6p-11f),
+ (float2) (0x1.1d4000p+0f, 0x1.0e62d0p-13f),
+ (float2) (0x1.204000p+0f, 0x1.1dc430p-11f),
+ (float2) (0x1.238000p+0f, 0x1.e9b9d4p-14f),
+ (float2) (0x1.268000p+0f, 0x1.a2b2f0p-11f),
+ (float2) (0x1.29c000p+0f, 0x1.4efa8ep-11f),
+ (float2) (0x1.2d0000p+0f, 0x1.42d372p-11f),
+ (float2) (0x1.304000p+0f, 0x1.7f0518p-11f),
+ (float2) (0x1.33c000p+0f, 0x1.164c82p-17f),
+ (float2) (0x1.370000p+0f, 0x1.a7373ap-12f),
+ (float2) (0x1.3a4000p+0f, 0x1.ed9a72p-11f),
+ (float2) (0x1.3dc000p+0f, 0x1.532608p-11f),
+ (float2) (0x1.414000p+0f, 0x1.0510fap-11f),
+ (float2) (0x1.44c000p+0f, 0x1.043030p-11f),
+ (float2) (0x1.484000p+0f, 0x1.515ae0p-11f),
+ (float2) (0x1.4bc000p+0f, 0x1.ed6a9ap-11f),
+ (float2) (0x1.4f8000p+0f, 0x1.b2769cp-12f),
+ (float2) (0x1.534000p+0f, 0x1.5ab4eap-15f),
+ (float2) (0x1.56c000p+0f, 0x1.a39b5ap-11f),
+ (float2) (0x1.5a8000p+0f, 0x1.83eea4p-11f),
+ (float2) (0x1.5e4000p+0f, 0x1.b78ad6p-11f),
+ (float2) (0x1.624000p+0f, 0x1.fac0e8p-14f),
+ (float2) (0x1.660000p+0f, 0x1.1c412ap-11f),
+ (float2) (0x1.6a0000p+0f, 0x1.3cccfep-13f),
+ (float2) (0x1.6dc000p+0f, 0x1.d91e32p-11f),
+ (float2) (0x1.71c000p+0f, 0x1.baf476p-11f),
+ (float2) (0x1.75c000p+0f, 0x1.f5ab20p-11f),
+ (float2) (0x1.7a0000p+0f, 0x1.1473eap-12f),
+ (float2) (0x1.7e0000p+0f, 0x1.799b66p-11f),
+ (float2) (0x1.824000p+0f, 0x1.89994cp-12f),
+ (float2) (0x1.868000p+0f, 0x1.b33688p-13f),
+ (float2) (0x1.8ac000p+0f, 0x1.ca8454p-13f),
+ (float2) (0x1.8f0000p+0f, 0x1.ae9914p-12f),
+ (float2) (0x1.934000p+0f, 0x1.9bd866p-11f),
+ (float2) (0x1.97c000p+0f, 0x1.829fdep-12f),
+ (float2) (0x1.9c4000p+0f, 0x1.230546p-13f),
+ (float2) (0x1.a0c000p+0f, 0x1.99ed76p-14f),
+ (float2) (0x1.a54000p+0f, 0x1.03b23ep-12f),
+ (float2) (0x1.a9c000p+0f, 0x1.35aabcp-11f),
+ (float2) (0x1.ae8000p+0f, 0x1.3f32b4p-13f),
+ (float2) (0x1.b30000p+0f, 0x1.d15c26p-11f),
+ (float2) (0x1.b7c000p+0f, 0x1.bb797cp-11f),
+ (float2) (0x1.bcc000p+0f, 0x1.e904bcp-16f),
+ (float2) (0x1.c18000p+0f, 0x1.9bdd84p-12f),
+ (float2) (0x1.c64000p+0f, 0x1.f8972ap-11f),
+ (float2) (0x1.cb4000p+0f, 0x1.906e76p-11f),
+ (float2) (0x1.d04000p+0f, 0x1.96a502p-11f),
+ (float2) (0x1.d58000p+0f, 0x1.8dcfbap-16f),
+ (float2) (0x1.da8000p+0f, 0x1.e603dap-12f),
+ (float2) (0x1.dfc000p+0f, 0x1.2e66f6p-13f),
+ (float2) (0x1.e50000p+0f, 0x1.773c58p-15f),
+ (float2) (0x1.ea4000p+0f, 0x1.5f4548p-13f),
+ (float2) (0x1.ef8000p+0f, 0x1.0df730p-11f),
+ (float2) (0x1.f50000p+0f, 0x1.d96db8p-14f),
+ (float2) (0x1.fa4000p+0f, 0x1.e0c0cep-11f),
+ (float2) (0x1.000000p+1f, 0x0.000000p+0f),
+};
+
TABLE_FUNCTION(float2, LOGE_TBL, loge_tbl);
TABLE_FUNCTION(float, LOG_INV_TBL, log_inv_tbl);
TABLE_FUNCTION(float2, LOG2_TBL, log2_tbl);
@@ -618,6 +754,8 @@ uint4 TABLE_MANGLE(pibits_tbl)(size_t idx) {
TABLE_FUNCTION(float2, SINHCOSH_TBL, sinhcosh_tbl);
TABLE_FUNCTION(float2, CBRT_TBL, cbrt_tbl);
+TABLE_FUNCTION(float, EXP_TBL, exp_tbl);
+TABLE_FUNCTION(float2, EXP_TBL_EP, exp_tbl_ep);
#ifdef cl_khr_fp64
diff --git a/generic/lib/math/tables.h b/generic/lib/math/tables.h
index c0b0ab1..1339a08 100644
--- a/generic/lib/math/tables.h
+++ b/generic/lib/math/tables.h
@@ -44,6 +44,8 @@ TABLE_FUNCTION_DECL(float2, log2_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);
+TABLE_FUNCTION_DECL(float2, exp_tbl_ep);
#ifdef cl_khr_fp64
--
2.9.3
More information about the Libclc-dev
mailing list