I recently read that some researchers were using the sigmoid-y Gudermannian function as an activation function, or at least investigating its utility. So I thought it would be cool to create an accurate and fast implementation of it for CUDA.
“The Guder… what?” I hear you say. It’s a slick little function, gd(x), that together with its inverse, gd-1(x), lets one move between circular and hyperbolic functions without going through complex numbers. While the mathematician Johann Heinrich Lambert (one of my mathematical heros, BTW) came up with a similar concept, he did not see any use for it and did not pursue it further. Christoph Gudermann (1798-1852) on the other hand fully developed the idea and published about it, because he thought it very useful for manipulating integrals and solving differential equations.
Gudermann did not name the function after himself. He was a humble man who worked as a high school teacher for nine years before becoming a university professor. Back then that was not so unusual, for example, Einstein initially got a teaching degree. The naming of the function occurred a decade after Gudermann’s death on the initiative of Arthur Cayley. Gudermann himself called these longitudinal functions: l for gd() and 𝔏 for gd-1(x).
Christoph Gudermann, “Theorie der Potential- oder cyklisch-hyperbolischen Functionen.” Part 2. Journal für die reine und angewandte Mathematik, Vol. 6, No. 2, 1830, pp. 162-194. (scan: https://babel.hathitrust.org/cgi/pt?id=mdp.39015026131113&seq=170)
There are various ways of computing the Gudermannian function from elementary functions, each with their own numerical quirks. One straightforward formula that follows directly from the definition is gd(x) = atan (sinh (x)).
As I thought about implementing it, I figured that the diminished exponential function expm1() might come in handy, and I thought I had the perfect way of implementing it efficiently with the help of the MUFU.EX2 and MUFU.LG2 instructions and Kahan’s famous trick [1]: Compute a preliminary result expm1(x) := exp(x) - 1, which obviously has a numerical issue with subtractive cancellation for x near 0, and then apply an adjustment factor if the function argument is close to the origin.
Kahan’s famous solution is to take the logarithm to compute the necessary adjustment. Sounds crazy on first encounter, but it works like a charm. People used this technique a lot before programming languages added a dedicated expm1() function to their standard math libraries. For C this happened with ISO-C99 and for C++ when ISO-C++11 incorporated the C99 math library by reference. Interesting little fact: I used this in the first release of CUDA – along with the similar trick for computing log1p() from log() – because I did not have enough time left for a proper implementation.
Massive numerical errors in my implementation soon reminded me that MUFU.LG2 is based on fixed-point computation and therefore has terrible relative error near 1, which is precisely where Kahan’s trick needs the logarithm to be accurate. I had to expend a fair number of additional instructions to supply my own little log2() implementation that uses a minimax polynomial approximation. In the end I wound up with my_hexpm1() below, that compiles to 19 machine instructions when using CUDA 13.3 to compile for an sm_89 platform. Only a single argument does not deliver a correctly rounded result in the regular path of computation and had to be fixed up as a special case at the end.
my_hexpm1(__half):
PRMT R4, R4, 0x7710, RZ // 1
MOV R3, 0x3e7a63c5 // 2
HADD2.F32 R4, -RZ, R4.H0_H0 // 3
FMUL R5, R4.reuse, 1.4426950216293334961 // 3
FSETP.GE.AND P0, PT, |R4|.reuse, 0.11400000005960464478, PT // 4
FSETP.GEU.AND P1, PT, |R4|, 0.00048828125, PT // 5
MUFU.EX2 R0, R5 // 6
FADD R0, R0, -1 // 7
FFMA R3, R0, -R3, 0.29260459542274475098 // 8
FFMA R3, R0, R3, -0.36064800620079040527 // 9
FFMA R3, R0, R3, 0.4808794558048248291 // 10
FFMA R3, R0, R3, -0.72134757041931152344 // 11
FFMA R3, R0, R3, 1.4426950216293334961 // 12
MUFU.RCP R6, R3 // 13
FMUL R7, R5, R6 // 14
@!P0 FSEL R0, R4.reuse, R7, !P1 // 15
FSETP.NEU.AND P0, PT, R4, -0.46435546875, PT // 16
F2F.F16.F32 R0, R0 // 17
SEL R4, R0, 0xb5f1, P0 // 19
RET.ABS.NODEC R20 0x0
My source code is as follows:
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
#include <cuda_fp16.h>
/*
Copyright 2026, Norbert Juffa
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.
*/
__forceinline__ __device__ float raw_ex2 (float a)
{
float r;
asm ("ex2.approx.ftz.f32 %0,%1;" : "=f"(r) : "f"(a));
return r;
}
__forceinline__ __device__ float raw_rcp (float a)
{
float r;
asm ("rcp.approx.ftz.f32 %0,%1;" : "=f"(r) : "f"(a));
return r;
}
__device__ __half my_hexpm1 (__half ha)
{
constexpr float L2E = 1.442695041f;
float aa, r, u, m, a;
__half hr;
aa = __half2float (ha);
a = __fmul_rn (aa, L2E);
u = raw_ex2 (a);
m = u - 1.0f; // preliminary result
// compute log2(1+m) = log2(u) for correction (alas, cannot use MUFU.LG2)
u = -2.445212113e-1f;
u = fmaf (u, m, 2.926046084e-1f);
u = fmaf (u, m, -3.606479971e-1f);
u = fmaf (u, m, 4.808794555e-1f);
u = fmaf (u, m, -7.213475602e-1f);
u = fmaf (u, m, 1.442695055e+0f);
r = raw_rcp (u) * a; // apply correction to preliminary result
if (fabsf (aa) < 0x1.0p-11f) r = aa; // expm1(x) = x near zero
if (fabsf (aa) >= 0.114f) r = m; // for larger |x| no correction is needed
hr = __float2half_rn (r);
// fix the one case that is not rounded correctly
if (aa==__half2float(__ushort_as_half(0xb76e))) hr=__ushort_as_half(0xb5f1);
return hr;
}
__global__ void test_hexpm1 (int *diffs)
{
int tid = blockDim.x * blockIdx.x + threadIdx.x;
unsigned short resi, refi, argi = tid;
__half arg, res, ref;
arg = __ushort_as_half (argi);
res = my_hexpm1 (arg);
ref = __double2half (expm1 ((double)__half2float (arg)));
resi = __half_as_ushort (res);
refi = __half_as_ushort (ref);
if (resi != refi) {
atomicAdd (diffs, 1);
printf ("hexpm1::arg=%04x (%11.4e) res=%04x (%11.4e) ref=%04x (%11.4e) %23.13a diff=%15.6a\n", argi, (float)arg, resi, (float)res, refi, (float)ref, expm1((double)__half2float (arg)), (float)res - (float)ref);
}
}
#define LEN (65536)
#define THREADS_PER_BLOCK (128)
int main (void)
{
dim3 dimBlock(THREADS_PER_BLOCK);
int threadBlocks = LEN / THREADS_PER_BLOCK;
dim3 dimGrid(threadBlocks);
int diffs = -1, *diffs_d = 0;
/* Select GPU to run on */
int dev = 0;
struct cudaDeviceProp props;
cudaSetDevice (dev);
cudaGetDeviceProperties (&props, dev);
printf ("test_hexpm1: running on device %d (%s)\n", dev, props.name);
cudaMalloc ((void**)&diffs_d, sizeof diffs);
cudaMemset (diffs_d, 0x00, sizeof diffs);
test_hexpm1<<<dimGrid,dimBlock>>>(diffs_d);
cudaMemcpy (&diffs, diffs_d, sizeof diffs, cudaMemcpyDeviceToHost);
printf ("diffs = %d\n", diffs);
cudaFree (diffs_d);
cudaDeviceSynchronize();
return EXIT_SUCCESS;
}
[1] W. Kahan, “Interval arithmetic options in the proposed IEEE floating point arithmetic standard”. In: Karl L. E. Nickel (ed.), Interval Mathematics 1980, New York: Academic Press 1980, pp. 99-128. A detailed numerical analysis of Kahan’s trick can be found in Nicholas J. Higham’s book Accuracy and Stability of Numerical Algorithms, Second Edition, SIAM 2002, in section 1.14.1 on pp. 19-21. An implementation of expm1() making use of Kahan’s technique might look as follows:
/*
The following is based on: W. Kahan, "Interval arithmetic options in the
proposed IEEE floating point arithmetic standard". In: Karl L. E. Nickel
(ed.), "Interval Mathematics 1980", New York: Academic Press 1980, pp.
99-128. See pages 110-111. For a detailed explanation see Nicholas J.
Higham, "Accuracy and Stability of Numerical Algorithms, Second Edition",
SIAM 2002. In section 1.14.1 on pages 19-21.
*/
double my_expm1 (double x)
{
double u, m;
u = exp (x);
m = u - 1.0;
if (m == 0.0) {
// x very close to zero
m = x;
} else if (fabs (x) < 1.0) {
// x somewhat close zero
u = log (u);
m = m * x;
m = m / u;
}
return m;
}