expm1(x) = exp(x) - 1, but accurate for |x| near zero. The function is most frequently used as a building block in financial computations and hyperbolic functions such as sinh().
When compared to CUDA 8’s built-in implementation, the implementation below achieves a throughput of 20.54e9 function evaluations per second on a Quadro K2200 (basically a GTX 750 Ti with larger memory) vs 18.85e9 for the implementation of CUDA 8, for a performance advantage of 9%. The maximum error is reduced from 1.44464 ulp in CUDA 8 to 1.00978 ulp. It would be nice if the function were faithfully rounded, i.e. maximum error < 1 ulp, but so far I have only managed to come close, with about a dozen inputs triggering error in excess of 1 ulp.
Given that expm1() isn’t all that frequently use, this is probably of limited utility, but after almost two years of tinkering with this every now and then I figured I might as well post the code. The code has been tested exhaustively for all possible inputs on a GPU with compute capability 5.0 using default compilation mode. Adjustments may be necessary when using other compilations modes or when porting to non-GPU platforms.
Code below superseded by the version in post #2 below.
/*
[s] Copyright (c) 2015-2017, Norbert Juffa
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
/* Compute exponential base e minus 1. Maximum ulp error = 1.00978
i = rint(a/log(2)), f = a-i*log(2). Then expm1(a) = 2**i * (expm1(f)+1) - 1.
Compute r = expm1(f). Then expm1(a)= 2 * (0.5 * 2**i * r + 0.5 * 2**i - 0.5).
t = 0.5*2**i. Thus expm1(a) = 2*(t + fma(r, t, -0.5)) or 2*fma(r, t, t-0.5).
When i == 1, expm1(a) = 2*(r + 0.5).
NOTE: Scale factor b is only applied if i < 0 or i > 1 (should be power of 2)
*/
__device__ float my_expm1f_scaled_unchecked (float a, float b)
{
float f, j, r, t, u;
int i;
// exp(a) = 2**i * exp(f); i = rintf (a / log(2))
j = fmaf (1.442695f, a, 12582912.f); // 0x1.715476p0, 0x1.8p23
i = __float_as_int (j);
j = j - 12582912.f; // 0x1.8p23
f = fmaf (j, -6.93145752e-1f, a); // -0x1.62e400p-1 // log_2_hi
f = fmaf (j, -1.42860677e-6f, f); // -0x1.7f7d1cp-20 // log_2_lo
// approximate r = exp(f)-1 on interval [-log(2)/2, +log(2)/2]
r = 1.98543072e-4f;
r = fmaf (r, f, 1.39400735e-3f);
r = fmaf (r, f, 8.33349116e-3f);
r = fmaf (r, f, 4.16664928e-2f);
r = fmaf (r, f, 1.66666716e-1f);
r = fmaf (r, f, 5.00000000e-1f);
t = (j == 1) ? (f + 0.5f) : f;
r = fmaf (r, f, 0.0f);
r = fmaf (r, f, t);
// Scale result from primary approximation range to full range
u = 0.5f * b;
t = __int_as_float ((i << 23) + __float_as_int (u)); // b*2**(i-1)
t = (j > 17) ? (fmaf (r, t, -u) + t) : fmaf (r, t, t - u);
if (j == 1) t = r;
if (j != 0) r = t + t;
return r;
}
/* Compute exponential base e minus 1. max ulp err = 1.00978 */
__device__ float my_expm1f (float a)
{
float r;
r = my_expm1f_scaled_unchecked (a, 1.0f);
/* handle severe overflow and underflow */
if (fabsf (a - 1.0f) > 88.0f) {
asm ("ex2.approx.ftz.f32 %0,%1;" : "=f"(r) : "f"(a));
r = fmaf (r, r, -1.0f);
}
return r;
}[/s]