As a follow-up to my recent weekend project, I checked whether the accuracy of CUDA’s built-in acosf() could be improved upon without negatively impacting performance. It turns out that this is possible, and one can even save one instruction in the process (although one might view this more as a compiler artifact than a proper micro-optimization).
CUDA’s current acosf() implementation compiles to 26 instructions for an sm_89 target, and achieves a maximum error < 1.34 ulps. my_acosf() below compiles to 25 instructions for an sm_89 target, and its maximum error is < 1.11 ulps.
my_acosf(float):
MOV R3, 0x3f000000
FSETP.GT.AND P1, PT, |R4|.reuse, 0.5625, PT
MOV R7, 0x3d41a000
FFMA R3, -|R4|, R3, 0.5
MUFU.RSQ R0, R3
FSETP.NEU.AND P0, PT, R3.reuse, RZ, PT
FMUL R6, R3, R0
FMUL R5, R0, 0.5
FFMA R0, R6, -R6, R3
@P0 FFMA R3, R5, R0, R6
FSETP.GTU.AND P0, PT, R4, 0.5625, PT
FSEL R5, R3, |R4|, P1
LOP3.LUT R5, R5, 0x80000000, R4, 0xf8, !PT
FSEL R5, R5, -R5, P1
FMUL R0, R5, R5
FFMA R3, R0, R7, 0.006988525390625
FFMA R3, R0, R3, 0.033894307911396026611
FFMA R3, R0, R3, 0.044259704649448394775
FFMA R3, R0, R3, 0.075017884373664855957
FFMA R3, R0, R3, 0.16666643321514129639
FMUL R0, R0, R3
@!P0 MOV R3, 0x3fd774eb
FFMA R4, R5, R0, R5
@!P0 FFMA R4, R3, 0.93318945169448852539, R4
@P1 FADD R4, R4, R4
RET.ABS.NODEC R20 0x0
In case someone finds this interesting, I am sharing my source code here:
/*
Copyright (c) 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_rsqrt (float a)
{
float r;
asm ("rsqrt.approx.ftz.f32 %0,%1;" : "=f"(r) : "f"(a));
return r;
}
__forceinline__ __device__ float copysignf_pos (float a, float b)
{
return __int_as_float (__float_as_int (a) | (__float_as_int (b) & 0x80000000));
}
// compute sqrt(a) from rsqrt approximation using Schoenhage's coupled iteration
__forceinline__ __device__ float sqrtf_schoenhage (float a)
{
float ra = raw_rsqrt (a);
float s0 = a * ra;
float r0 = 0.5f * ra;
float s1 = fmaf (fmaf (s0, -s0, a), r0, s0);
return (a == 0) ? a : s1;
}
// compute acos(x) with maximum error of < 1.11 ulp
__device__ float my_acosf (float x)
{
const float SWITCHOVER = 0.5625f;
const float pio2_fac1 = 0.93318945f;
const float pio2_fac2 = 1.68325555f;
float res, sq, redx, ax = fabsf (x);
// arccos(x) = 2 * arcsin (sqrt((1 - x) / 2))
redx = sqrtf_schoenhage (fmaf (0.5f, -ax, 0.5f));
redx = (ax > SWITCHOVER) ? redx : ax;
redx = copysignf_pos (redx, x);
redx = (ax > SWITCHOVER) ? redx : (-redx);
// approximate asin(redx) on [0, SWITCHOVER]
sq = redx * redx;
res = 0x1.834000p-5f;
res = fmaf (res, sq, 0x1.ca0000p-8f);
res = fmaf (res, sq, 0x1.15a984p-5f);
res = fmaf (res, sq, 0x1.6a9354p-5f);
res = fmaf (res, sq, 0x1.3345f4p-4f);
res = fmaf (res, sq, 0x1.555536p-3f);
res = res * sq;
res = fmaf (res, redx, redx);
// map asin(redx) to acos(x) according to x and magnitude of |x|
res = (x <= SWITCHOVER) ? fmaf (pio2_fac1, pio2_fac2, res) : res;
res = (ax > SWITCHOVER) ? (res + res) : res;
return res;
}