275 lines · plain
1//===----------------------------------------------------------------------===//2//3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.4// See https://llvm.org/LICENSE.txt for license information.5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception6//7//===----------------------------------------------------------------------===//8 9// This version is derived from the generic fma software implementation10// (__clc_sw_fma), but avoids the use of ulong in favor of uint2. The logic has11// been updated as appropriate.12 13#include <clc/clc_as_type.h>14#include <clc/float/definitions.h>15#include <clc/integer/clc_abs.h>16#include <clc/integer/clc_clz.h>17#include <clc/integer/clc_hadd.h>18#include <clc/integer/clc_mul_hi.h>19#include <clc/integer/definitions.h>20#include <clc/math/clc_mad.h>21#include <clc/math/math.h>22#include <clc/relational/clc_isinf.h>23#include <clc/relational/clc_isnan.h>24#include <clc/shared/clc_max.h>25 26struct fp {27 uint2 mantissa;28 int exponent;29 uint sign;30};31 32static uint2 u2_set(uint hi, uint lo) {33 uint2 res;34 res.lo = lo;35 res.hi = hi;36 return res;37}38 39static uint2 u2_set_u(uint val) { return u2_set(0, val); }40 41static uint2 u2_mul(uint a, uint b) {42 uint2 res;43 res.hi = __clc_mul_hi(a, b);44 res.lo = a * b;45 return res;46}47 48static uint2 u2_sll(uint2 val, uint shift) {49 if (shift == 0)50 return val;51 if (shift < 32) {52 val.hi <<= shift;53 val.hi |= val.lo >> (32 - shift);54 val.lo <<= shift;55 } else {56 val.hi = val.lo << (shift - 32);57 val.lo = 0;58 }59 return val;60}61 62static uint2 u2_srl(uint2 val, uint shift) {63 if (shift == 0)64 return val;65 if (shift < 32) {66 val.lo >>= shift;67 val.lo |= val.hi << (32 - shift);68 val.hi >>= shift;69 } else {70 val.lo = val.hi >> (shift - 32);71 val.hi = 0;72 }73 return val;74}75 76static uint2 u2_or(uint2 a, uint b) {77 a.lo |= b;78 return a;79}80 81static uint2 u2_and(uint2 a, uint2 b) {82 a.lo &= b.lo;83 a.hi &= b.hi;84 return a;85}86 87static uint2 u2_add(uint2 a, uint2 b) {88 uint carry = (__clc_hadd(a.lo, b.lo) >> 31) & 0x1;89 a.lo += b.lo;90 a.hi += b.hi + carry;91 return a;92}93 94static uint2 u2_add_u(uint2 a, uint b) { return u2_add(a, u2_set_u(b)); }95 96static uint2 u2_inv(uint2 a) {97 a.lo = ~a.lo;98 a.hi = ~a.hi;99 return u2_add_u(a, 1);100}101 102static uint u2_clz(uint2 a) {103 uint leading_zeroes = __clc_clz(a.hi);104 if (leading_zeroes == 32) {105 leading_zeroes += __clc_clz(a.lo);106 }107 return leading_zeroes;108}109 110static bool u2_eq(uint2 a, uint2 b) { return a.lo == b.lo && a.hi == b.hi; }111 112static bool u2_zero(uint2 a) { return u2_eq(a, u2_set_u(0)); }113 114static bool u2_gt(uint2 a, uint2 b) {115 return a.hi > b.hi || (a.hi == b.hi && a.lo > b.lo);116}117 118_CLC_DEF _CLC_OVERLOAD float __clc_sw_fma(float a, float b, float c) {119 /* special cases */120 if (__clc_isnan(a) || __clc_isnan(b) || __clc_isnan(c) || __clc_isinf(a) ||121 __clc_isinf(b)) {122 return __clc_mad(a, b, c);123 }124 125 /* If only c is inf, and both a,b are regular numbers, the result is c*/126 if (__clc_isinf(c)) {127 return c;128 }129 130 a = __clc_flush_denormal_if_not_supported(a);131 b = __clc_flush_denormal_if_not_supported(b);132 c = __clc_flush_denormal_if_not_supported(c);133 134 if (a == 0.0f || b == 0.0f) {135 return c;136 }137 138 if (c == 0) {139 return a * b;140 }141 142 struct fp st_a, st_b, st_c;143 144 st_a.exponent = a == .0f ? 0 : ((__clc_as_uint(a) & 0x7f800000) >> 23) - 127;145 st_b.exponent = b == .0f ? 0 : ((__clc_as_uint(b) & 0x7f800000) >> 23) - 127;146 st_c.exponent = c == .0f ? 0 : ((__clc_as_uint(c) & 0x7f800000) >> 23) - 127;147 148 st_a.mantissa =149 u2_set_u(a == .0f ? 0 : (__clc_as_uint(a) & 0x7fffff) | 0x800000);150 st_b.mantissa =151 u2_set_u(b == .0f ? 0 : (__clc_as_uint(b) & 0x7fffff) | 0x800000);152 st_c.mantissa =153 u2_set_u(c == .0f ? 0 : (__clc_as_uint(c) & 0x7fffff) | 0x800000);154 155 st_a.sign = __clc_as_uint(a) & 0x80000000;156 st_b.sign = __clc_as_uint(b) & 0x80000000;157 st_c.sign = __clc_as_uint(c) & 0x80000000;158 159 // Multiplication.160 // Move the product to the highest bits to maximize precision161 // mantissa is 24 bits => product is 48 bits, 2bits non-fraction.162 // Add one bit for future addition overflow,163 // add another bit to detect subtraction underflow164 struct fp st_mul;165 st_mul.sign = st_a.sign ^ st_b.sign;166 st_mul.mantissa = u2_sll(u2_mul(st_a.mantissa.lo, st_b.mantissa.lo), 14);167 st_mul.exponent =168 !u2_zero(st_mul.mantissa) ? st_a.exponent + st_b.exponent : 0;169 170 // FIXME: Detecting a == 0 || b == 0 above crashed GCN isel171 if (st_mul.exponent == 0 && u2_zero(st_mul.mantissa))172 return c;173 174// Mantissa is 23 fractional bits, shift it the same way as product mantissa175#define C_ADJUST 37ul176 177 // both exponents are bias adjusted178 int exp_diff = st_mul.exponent - st_c.exponent;179 180 st_c.mantissa = u2_sll(st_c.mantissa, C_ADJUST);181 uint2 cutoff_bits = u2_set_u(0);182 uint2 cutoff_mask = u2_add(u2_sll(u2_set_u(1), __clc_abs(exp_diff)),183 u2_set(0xffffffff, 0xffffffff));184 if (exp_diff > 0) {185 cutoff_bits =186 exp_diff >= 64 ? st_c.mantissa : u2_and(st_c.mantissa, cutoff_mask);187 st_c.mantissa =188 exp_diff >= 64 ? u2_set_u(0) : u2_srl(st_c.mantissa, exp_diff);189 } else {190 cutoff_bits = -exp_diff >= 64 ? st_mul.mantissa191 : u2_and(st_mul.mantissa, cutoff_mask);192 st_mul.mantissa =193 -exp_diff >= 64 ? u2_set_u(0) : u2_srl(st_mul.mantissa, -exp_diff);194 }195 196 struct fp st_fma;197 st_fma.sign = st_mul.sign;198 st_fma.exponent = __clc_max(st_mul.exponent, st_c.exponent);199 if (st_c.sign == st_mul.sign) {200 st_fma.mantissa = u2_add(st_mul.mantissa, st_c.mantissa);201 } else {202 // cutoff bits borrow one203 st_fma.mantissa =204 u2_add(u2_add(st_mul.mantissa, u2_inv(st_c.mantissa)),205 (!u2_zero(cutoff_bits) && (st_mul.exponent > st_c.exponent)206 ? u2_set(0xffffffff, 0xffffffff)207 : u2_set_u(0)));208 }209 210 // underflow: st_c.sign != st_mul.sign, and magnitude switches the sign211 if (u2_gt(st_fma.mantissa, u2_set(0x7fffffff, 0xffffffff))) {212 st_fma.mantissa = u2_inv(st_fma.mantissa);213 st_fma.sign = st_mul.sign ^ 0x80000000;214 }215 216 // detect overflow/underflow217 int overflow_bits = 3 - u2_clz(st_fma.mantissa);218 219 // adjust exponent220 st_fma.exponent += overflow_bits;221 222 // handle underflow223 if (overflow_bits < 0) {224 st_fma.mantissa = u2_sll(st_fma.mantissa, -overflow_bits);225 overflow_bits = 0;226 }227 228 // rounding229 uint2 trunc_mask = u2_add(u2_sll(u2_set_u(1), C_ADJUST + overflow_bits),230 u2_set(0xffffffff, 0xffffffff));231 uint2 trunc_bits =232 u2_or(u2_and(st_fma.mantissa, trunc_mask), !u2_zero(cutoff_bits));233 uint2 last_bit =234 u2_and(st_fma.mantissa, u2_sll(u2_set_u(1), C_ADJUST + overflow_bits));235 uint2 grs_bits = u2_sll(u2_set_u(4), C_ADJUST - 3 + overflow_bits);236 237 // round to nearest even238 if (u2_gt(trunc_bits, grs_bits) ||239 (u2_eq(trunc_bits, grs_bits) && !u2_zero(last_bit))) {240 st_fma.mantissa =241 u2_add(st_fma.mantissa, u2_sll(u2_set_u(1), C_ADJUST + overflow_bits));242 }243 244 // Shift mantissa back to bit 23245 st_fma.mantissa = u2_srl(st_fma.mantissa, C_ADJUST + overflow_bits);246 247 // Detect rounding overflow248 if (u2_gt(st_fma.mantissa, u2_set_u(0xffffff))) {249 ++st_fma.exponent;250 st_fma.mantissa = u2_srl(st_fma.mantissa, 1);251 }252 253 if (u2_zero(st_fma.mantissa)) {254 return 0.0f;255 }256 257 // Flating point range limit258 if (st_fma.exponent > 127) {259 return __clc_as_float(__clc_as_uint(INFINITY) | st_fma.sign);260 }261 262 // Flush denormals263 if (st_fma.exponent <= -127) {264 return __clc_as_float(st_fma.sign);265 }266 267 return __clc_as_float(st_fma.sign | ((st_fma.exponent + 127) << 23) |268 ((uint)st_fma.mantissa.lo & 0x7fffff));269}270 271#define __CLC_FLOAT_ONLY272#define __CLC_FUNCTION __clc_sw_fma273#define __CLC_BODY <clc/shared/ternary_def_scalarize.inc>274#include <clc/math/gentype.inc>275