I noticed that CUDA 13.4 does not offer any half-precision inverse trigonometric functions. Since the arctangent is used as an activation function, that seemed like the one CUDA programmers would most likely need.
The development process here followed my (by now) usual path. I grabbed working FP32 code of mine, slightly simplified it for the lower accuracy requirement, then tried to minimize the number of arguments that do not result in a correctly rounded FP16 result by the judicious use of directed rounded modes. I got down to one argument needing manual fixup.
Compiled for an sm_89 target, my_hatan() compiles to 20 instructions. I looked at some implementation alternatives but could not any that required fewer machine instructions. Maybe someone else has a cool idea how to compute this more efficiently.
#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.
*/
/* Use MUFU.RCP directly */
__forceinline__ __device__ float gpu_rcp_approx (float divisor)
{
float r;
asm ("rcp.approx.ftz.f32 %0,%1;\n\t" : "=f"(r) : "f"(divisor));
return r;
}
/* Use MUFU.RSQ directly */
__forceinline__ __device__ float gpu_rsqrt_approx (float divisor)
{
float r;
asm ("rsqrt.approx.ftz.f32 %0,%1;\n\t" : "=f"(r) : "f"(divisor));
return r;
}
/* Use MUFU.SQRT directly */
__forceinline__ __device__ float gpu_sqrt_approx (float radicand)
{
float r;
asm ("sqrt.approx.ftz.f32 %0,%1;\n\t" : "=f"(r) : "f"(radicand));
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__ float atan_poly (float a)
{
float r, s;
s = a * a;
r = -0x1.3c0cb6p-8f ;
r = fmaf (r, s, 0x1.953e30p-6f);
r = fmaf (r, s, -0x1.ed2f1ep-5f);
r = fmaf (r, s, 0x1.984f04p-4f);
r = fmaf (r, s, -0x1.1f9106p-3f);
r = fmaf (r, s, 0x1.991260p-3f);
r = fmaf (r, s, -0x1.5552e8p-2f);
r = r * s;
r = __fmaf_ru (r, a, a);
return r;
}
__device__ float atanf_core (float a)
{
float r, t;
t = fabsf (a);
r = t;
if (t > 1.0f) r = gpu_rcp_approx (r);
r = atan_poly (r);
if (t > 1.0f) r = __fadd_rz (1.57079633f, -r);
// one argument |a| does not produce a correctly rounded result; fix it here
if (t == __half2float(__ushort_as_half(0x466e))) {
r = __half2float(__ushort_as_half(0x3dab));
}
r = my_copysignf (r, a);
return r;
}
__device__ __half my_hatan (__half ha)
{
__half hr;
float r, a = __half2float (ha);
r = atanf_core (a);
hr = __float2half_rn (r);
return hr;
}
__global__ void test_hatan (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_hatan (arg);
ref = __double2half (atan((double)__half2float (arg)));
resi = __half_as_ushort (res);
refi = __half_as_ushort (ref);
if (resi != refi) {
atomicAdd (diffs, 1);
printf ("hatan::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, atan((double)__half2float (arg)), (float)res-(float)ref);
}
}
// Macro to catch CUDA errors in CUDA runtime calls
#define CUDA_SAFE_CALL(call) \
do { \
cudaError_t err = call; \
if (cudaSuccess != err) { \
fprintf (stderr, "Cuda error in file '%s' in line %i : %s.\n",\
__FILE__, __LINE__, cudaGetErrorString(err) ); \
exit(EXIT_FAILURE); \
} \
} while (0)
// Macro to catch CUDA errors in kernel launches
#define CHECK_LAUNCH_ERROR() \
do { \
/* Check synchronous errors, i.e. pre-launch */ \
cudaError_t err = cudaGetLastError(); \
if (cudaSuccess != err) { \
fprintf (stderr, "Cuda error in file '%s' in line %i : %s.\n",\
__FILE__, __LINE__, cudaGetErrorString(err) ); \
exit(EXIT_FAILURE); \
} \
/* Check asynchronous errors, i.e. kernel failed (ULF) */ \
err = cudaDeviceSynchronize(); \
if (cudaSuccess != err) { \
fprintf (stderr, "Cuda error in file '%s' in line %i : %s.\n",\
__FILE__, __LINE__, cudaGetErrorString( err) ); \
exit(EXIT_FAILURE); \
} \
} while (0)
#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, *diffs_d = 0;
/* Select GPU to run on */
int dev = 0;
struct cudaDeviceProp props;
CUDA_SAFE_CALL (cudaSetDevice (dev));
CUDA_SAFE_CALL (cudaGetDeviceProperties (&props, dev));
printf ("test_hatan: running on device %d (%s)\n", dev, props.name);
CUDA_SAFE_CALL (cudaMalloc ((void**)&diffs_d, sizeof diffs));
CUDA_SAFE_CALL (cudaMemset (diffs_d, 0x00, sizeof diffs));
test_hatan<<<dimGrid,dimBlock>>>(diffs_d);
CHECK_LAUNCH_ERROR();
CUDA_SAFE_CALL (cudaMemcpy (&diffs, diffs_d, sizeof diffs,
cudaMemcpyDeviceToHost));
printf ("diffs = %d\n", diffs);
CUDA_SAFE_CALL (cudaFree (diffs_d));
cudaDeviceSynchronize();
return EXIT_SUCCESS;
}