For some exploratory work I was looking for a highly efficient implementation of arc cosine with an accuracy of 16 bits that is significantly faster than CUDA’s built-in acosf() implementation. For reference, acosf() from the CUDA 12.8 standard math library achieves a maximum error of 1.34 ulps (FP32) and compiles to 26 instructions for an sm_89 target.
Given my positive experience with using Borchardt’s algorithm (BA) with convergence acceleration for computing arc tangent with reduced precision, I tried that first. I found the resulting code suitable for my purposes, so I am sharing it here. I did not quite get to 16-bit accuracy but close enough. Using Borchardt’s algorithm, we have acos(x) = √(1-x2)/BA(x,1). I used a heuristic optimizer to find “optimal” coefficients for the convergence acceleration. The resulting implementation acosf_borchardt_15() compiles to 12 instructions for an sm_89 target. The direct use of the MUFU.SQRT instruction limits the code to compute capability 5.2 (second generation Maxwell) or higher. Given that CC 5.2 devices first shipped in 2014, this should not be a significant limitation. acosf_borchardt_15() limits maximum relative error to 1.9514e-5, and its maximum error is 303.08 ulps (FP32).
acosf_borchardt_15(float):
MOV R3, 0x3f000000
FFMA R0, R4.reuse, -R4, 1
MOV R8, 0x3ea39000
FSETP.GEU.AND P0, PT, R4.reuse, RZ, PT
FFMA R3, |R4|, R3, 0.5
MUFU.SQRT R0, R0
FFMA R5, R3, R8, -0.0148525238037109375
MUFU.SQRT R6, R3
FFMA R5, R6, 0.69541257619857788086, R5
MUFU.RCP R5, R5
@P0 FMUL R4, R0.reuse, R5.reuse
@!P0 FFMA R4, -R0, R5, 3.1415927410125732422
RET.ABS.NODEC R20 0x0
While I have no immediate use for even lower accuracy implementations of arc cosine, other forums participants might be interested in that. I tried Borchardt’s algorithm with no iterations but convergence acceleration, again using a heuristic optimizer to puzzle out “optimal” coefficients. The resulting implementation, acosf_borchardt_7() requires just eight instructions on sm_89 and is accurate to about 7.5 bits. Specifically, its maximum relative error is 5.4812e-3, while the maximum error is 91951 ulps (FP32).
acosf_borchardt_7(float):
MOV R3, 0x3ebb136f
FFMA R0, R4.reuse, -R4, 1
FSETP.GEU.AND P0, PT, R4, RZ, PT
FFMA R3, |R4|, R3, 0.64012819528579711914
MUFU.SQRT R0, R0
MUFU.RCP R3, R3
@P0 FMUL R4, R0.reuse, R3.reuse
@!P0 FFMA R4, -R0, R3, 3.1415927410125732422
RET.ABS.NODEC R20 0x0
An alternative implementation based on Robert Harley’s observation that acos(x) ≈ √(2*(1-x)) combined with a linear correction achieves similar accuracy and also compiles to eight instructions for an sm_89 target. This might be preferable in various contexts due to lower MUFU utilization. The maximum relative error is 5.1408e-3, while the maximum error is 85791 ulps (FP32).
acosf_harley(float):
MOV R3, 0x40000000
FSETP.GEU.AND P0, PT, R4, RZ, PT
FFMA R0, -|R4|, R3, 2
MOV R3, 0x3dd70ff6
MUFU.SQRT R5, R0
FFMA R4, -|R4|, R3, 0.10501091182231903076
FFMA R4, R5, R4, R5
@!P0 FADD R4, -R4, 3.1415927410125732422
RET.ABS.NODEC R20 0x0
/*
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_rcp (float a)
{
asm ("rcp.approx.ftz.f32 %0,%0;" : "+f"(a));
return a;
}
__forceinline__ __device__ float raw_sqrt (float a)
{
asm ("sqrt.approx.ftz.f32 %0,%0;" : "+f"(a));
return a;
}
/* Compute acos(x) using Borchardt's algorithm with convergence acceleration.
Approximation is accurate to about 15.5 bits. The maximum relative error
is 1.9514e-5, while the maximum error is 303.08 ulp.
*/
__device__ float acosf_borchardt_15 (float x)
{
const float PI = 3.14159265f;
float a, b, r, s, t;
t = fmaf (x, -x, 1.0f);
a = fmaf (0.5f, fabsf (x), 0.5f);
b = raw_sqrt (a);
r = fmaf (3.19458008e-1f, a, -1.48525238e-2f);
s = raw_sqrt (t);
r = fmaf (6.95412576e-1f, b, r);
r = raw_rcp (r);
r = (x < 0) ? fmaf (-s, r, PI) : (s * r);
return r;
}
/* Compute acos(x) using Borchardt's algorithm with convergence acceleration.
Approximation is accurate to about 7.5 bits. The maximum relative error is
5.4812e-3, while the maximum error is 91951 ulp.
*/
__device__ float acosf_borchardt_7 (float x)
{
const float PI = 3.14159265f;
float r, s, t;
t = fmaf (x, -x, 1.0f);
s = raw_sqrt (t);
r = fmaf (3.65382433e-1f, fabsf (x), 6.40128374e-1f);
r = raw_rcp (r);
r = (x < 0) ? fmaf (-s, r, PI) : (s * r);
return r;
}
/*
Compute acos(x) using an observation by Robert Harley in comp.arch.arithmetic
on 1996/07/12.
https://groups.google.com/forum/#!original/comp.arch.arithmetic/wqCPkCCXqWs/T9qCkHtGE2YJ.
Approximation is accurate to about 7.5 bits. The maximum relative error is
5.1408e-3, while the maximum error is 85791 ulp.
*/
__device__ float acosf_harley (float x)
{
const float PI = 3.14159265f;
const float C = 0.10501085f;
float r, s, t;
t = -fabsf (x);
s = raw_sqrt (fmaf (2, t, 2));
r = fmaf (fmaf (C, t, C), s, s);
if (x < 0) r = PI - r;
return r;
}