Since the hyperbolic tangent is a fairly popular activation function, I created a correctly rounded half-precision implementation for it from some existing code I had sitting around. This turned out to produce 18 machine instructions for an sm_89 target, compared to 20 machine instructions for the built-in `htanh() of CUDA 13.3, for a saving of two instructions.
Briefly looking over the SASS of the CUDA implementation, it seems to be implementing roughly the same algorithm, so these savings likely come down to implementation details. I did not bother to investigate further.
Compiling with CUDA 13.3 for an sm_89 target results in the following SASS for my_htanh():
my_htanh(__half):
PRMT R4, R4, 0x7710, RZ // 1
MOV R0, 0x3c835000 // 2
MOV R6, 0x40000000 // 3
HADD2.F32 R3, -RZ, R4.H0_H0 // 4
FMUL R4, |R3|.reuse, 2.8853900432586669922 // 5
FMUL R7, R3.reuse, R3 // 6
FSETP.GT.AND P0, PT, |R3|, 0.57421875, PT // 7
FFMA R0, R7.reuse, R0, -0.052451841533184051514 // 8
MUFU.EX2 R4, R4 // 9
FFMA R0, R7, R0, 0.13317251205444335938 // 10
FFMA R0, R7, R0, -0.33332800865173339844 // 11
FFMA R0, R7, R0, RZ // 12
FADD R5, R4, 1 // 13
MUFU.RCP R5, R5 // 14
FFMA R6, R5, -R6, 1 // 15
LOP3.LUT R4, R6, 0x80000000, R3.reuse, 0xb8, !PT // 16
@!P0 FFMA R4, R3, R0, R3 // 17
F2F.F16.F32 R4, R4 // 18
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;
}
// (a & ~b) | (c & b)
__forceinline__ __device__ unsigned int lop3_b8 (unsigned int a, unsigned int b, unsigned int c)
{
unsigned int r;
asm ("lop3.b32 %0,%1,%2,%3,0xb8;\n\t" : "=r"(r) : "r"(a), "r"(b), "r"(c));
return r;
}
__forceinline__ __device__ float my_copysignf (float x, float y)
{
unsigned int a = __float_as_int (x);
unsigned int c = __float_as_int (y);
unsigned int b = 0x80000000;
unsigned int r = lop3_b8 (a, b, c);
return __int_as_float (r);
}
__device__ __half my_htanh (__half ha)
{
float a, r, s, t, u;
__half hr;
a = __half2float (ha);
t = fabsf (a);
u = raw_rcp (1.0f + raw_ex2 (2.885390082f * t)); // 2*log2(e)
u = fmaf (-2.0f, u, 1.0f);
u = my_copysignf (u, a);
s = a * a;
r = 0x1.06a000p-6f;
r = fmaf (r, s, -0x1.adaf7cp-5f);
r = fmaf (r, s, 0x1.10bcc0p-3f);
r = fmaf (r, s, -0x1.5553f0p-2f);
r = fmaf (r, s, 0.0f);
r = fmaf (r, a, a);
if (t > 147.0f/256.0f) r = u;
hr = __float2half (r);
return hr;
}
__device__ __half htanh_wrapper (__half a)
{
return htanh (a);
}
__global__ void test_htanh (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_htanh (arg);
ref = __double2half (tanh((double)__half2float (arg)));
resi = __half_as_ushort (res);
refi = __half_as_ushort (ref);
if (resi != refi) {
atomicAdd (diffs, 1);
printf ("htanh::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, tanh((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_htanh: running on device %d (%s)\n", dev, props.name);
cudaMalloc ((void**)&diffs_d, sizeof diffs);
cudaMemset (diffs_d, 0x00, sizeof diffs);
test_htanh<<<dimGrid,dimBlock>>>(diffs_d);
cudaMemcpy (&diffs, diffs_d, sizeof diffs, cudaMemcpyDeviceToHost);
printf ("diffs = %d\n", diffs);
cudaFree (diffs_d);
cudaDeviceSynchronize();
return EXIT_SUCCESS;
}