I was toying with the idea how to use the SFU instructions for half-precision math functions of good accuracy (with a 11-bit mantissa, you can’t afford many ulps of error). By performing all intermediate computation in single-precision, one can side-step issues with underflow and overflow in intermediate computation, and the known numerical issues with use of the SFU instructions seem to mostly disappear once the final result is converted to half precision. The SFU instructions also provide decent special case handling, which saves instructions for handling those.
Below is one example, a half-precision version of sincos(). Note that the __fp16 type used is my own ‘typedef’; use the half-precision type name of your choice there. One nasty surprise was that there seems to be a consistent bias in either the argument reduction instruction RRO.SINCOS, or the MUFU.SIN and MUFU.COS instructions themselves (maybe due to an internal computation that is truncating instead of rounding?). So I had to add a bias-adjustment to the reduced argument for accurate results, which unfortunately currently adds three instructions.
Two variants are shown, with the more accurate one providing results that differ by at most 2 ulps from the correctly rounded results, while the less accurate one provides results that differ by at most 4 ulps from correctly rounded results, the worst case errors for the latter occur at:
[sin] error: arg= 1.06500000e+3 6429 res=-9.06586647e-5 85f1 ref=-9.04330600e-5 85ed
[sin] error: arg=-1.06500000e+3 e429 res= 9.06586647e-5 05f1 ref= 9.04330600e-5 05ed
/*
Copyright (c) 2016, Norbert Juffa
All rights reserved.
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_sin (float a)
{
float r;
asm ("sin.approx.ftz.f32 %0,%1;" : "=f"(r) : "f"(a));
return r;
}
__forceinline__ __device__ float raw_cos (float a)
{
float r;
asm ("cos.approx.ftz.f32 %0,%1;" : "=f"(r) : "f"(a));
return r;
}
__forceinline__ __device__ float copysignf_pos (float a, float b)
{
float r;
r = __int_as_float (__float_as_int (a) | (__float_as_int (b) & 0x80000000));
return r;
}
/* Use 'sf' suffix as per proposal N2016 in ISO/IEC JTC1 SC22 WG14 */
__device__ void my_sincossf (__fp16 ah, __fp16 *sp, __fp16 *cp)
{
float a, i, j, r, s, c;
const float BIAS_ADJ = 8.8e-8f;
a = __half2float (ah);
#if MORE_ACCURATE // maximal difference from correctly rounded result: 2 ulps
i = fmaf (a, 0.636619747f/2, 12582912.0f); // 1/pi, 0x1.8p+23
j = i - 12582912.0f; // 0x1.8p+23
r = fmaf (j, -3.14159203e+0f, a); // -0x1.921fb0p+01 // pi_high
r = fmaf (j, -6.27832947e-7f, r); // -0x1.5110b4p-21 // pi_low
r = r + copysignf_pos (BIAS_ADJ, r); // correct bias of trig func intrinsics
c = raw_cos (r);
s = raw_sin (r);
if (__float_as_int (i) & 1) {
s = 0.0f - s; // don't change "sign" of NaNs or create negative zeros
c = 0.0f - c; // don't change "sign" of NaNs or create negative zeros
}
#else // maximal difference from correctly rounded result: 3 ulps
i = fmaf (a, 0.636619747f/4, 12582912.0f); // 1/(2*pi), 0x1.8p+23
j = i - 12582912.0f; // 0x1.8p+23
r = fmaf (j, -6.28318405e+0f, a); // -0x1.921fb0p+02 // 2pi_high
r = fmaf (j, -1.25566589e-6f, r); // -0x1.5110b4p-20 // 2pi_low
r = r + copysignf_pos (BIAS_ADJ, r); // correct bias of trig func intrinsics
c = raw_cos (r);
s = raw_sin (r);
#endif
*sp = __float2half_rn (s);
*cp = __float2half_rn (c);
}