1/*
2 * Copyright (c) 2014 Advanced Micro Devices, Inc.
3 *
4 * Permission is hereby granted, free of charge, to any person obtaining a copy
5 * of this software and associated documentation files (the "Software"), to deal
6 * in the Software without restriction, including without limitation the rights
7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8 * copies of the Software, and to permit persons to whom the Software is
9 * furnished to do so, subject to the following conditions:
10 *
11 * The above copyright notice and this permission notice shall be included in
12 * all copies or substantial portions of the Software.
13 *
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
20 * THE SOFTWARE.
21 */
22
23// This version is derived from the generic fma software implementation
24// (__clc_sw_fma), but avoids the use of ulong in favor of uint2. The logic has
25// been updated as appropriate.
26
27#include <clc/clc.h>
28#include "../../../generic/lib/clcmacro.h"
29#include "../../../generic/lib/math/math.h"
30
31struct fp {
32  uint2 mantissa;
33  int exponent;
34  uint sign;
35};
36
37_CLC_DEF _CLC_OVERLOAD float fma(float a, float b, float c) {
38  /* special cases */
39  if (isnan(a) || isnan(b) || isnan(c) || isinf(a) || isinf(b)) {
40    return mad(a, b, c);
41  }
42
43  /* If only c is inf, and both a,b are regular numbers, the result is c*/
44  if (isinf(c)) {
45    return c;
46  }
47
48  a = __clc_flush_denormal_if_not_supported(a);
49  b = __clc_flush_denormal_if_not_supported(b);
50  c = __clc_flush_denormal_if_not_supported(c);
51
52  if (a == 0.0f || b == 0.0f) {
53    return c;
54  }
55
56  if (c == 0) {
57    return a * b;
58  }
59
60  struct fp st_a, st_b, st_c;
61
62  st_a.exponent = a == .0f ? 0 : ((as_uint(a) & 0x7f800000) >> 23) - 127;
63  st_b.exponent = b == .0f ? 0 : ((as_uint(b) & 0x7f800000) >> 23) - 127;
64  st_c.exponent = c == .0f ? 0 : ((as_uint(c) & 0x7f800000) >> 23) - 127;
65
66  st_a.mantissa.lo = a == .0f ? 0 : (as_uint(a) & 0x7fffff) | 0x800000;
67  st_b.mantissa.lo = b == .0f ? 0 : (as_uint(b) & 0x7fffff) | 0x800000;
68  st_c.mantissa.lo = c == .0f ? 0 : (as_uint(c) & 0x7fffff) | 0x800000;
69  st_a.mantissa.hi = 0;
70  st_b.mantissa.hi = 0;
71  st_c.mantissa.hi = 0;
72
73  st_a.sign = as_uint(a) & 0x80000000;
74  st_b.sign = as_uint(b) & 0x80000000;
75  st_c.sign = as_uint(c) & 0x80000000;
76
77  // Multiplication.
78  // Move the product to the highest bits to maximize precision
79  // mantissa is 24 bits => product is 48 bits, 2bits non-fraction.
80  // Add one bit for future addition overflow,
81  // add another bit to detect subtraction underflow
82  struct fp st_mul;
83  st_mul.sign = st_a.sign ^ st_b.sign;
84  st_mul.mantissa.hi = mul_hi(st_a.mantissa.lo, st_b.mantissa.lo);
85  st_mul.mantissa.lo = st_a.mantissa.lo * st_b.mantissa.lo;
86  uint upper_14bits = (st_mul.mantissa.lo >> 18) & 0x3fff;
87  st_mul.mantissa.lo <<= 14;
88  st_mul.mantissa.hi <<= 14;
89  st_mul.mantissa.hi |= upper_14bits;
90  st_mul.exponent = (st_mul.mantissa.lo != 0 || st_mul.mantissa.hi != 0)
91                        ? st_a.exponent + st_b.exponent
92                        : 0;
93
94// Mantissa is 23 fractional bits, shift it the same way as product mantissa
95#define C_ADJUST 37ul
96
97  // both exponents are bias adjusted
98  int exp_diff = st_mul.exponent - st_c.exponent;
99
100  uint abs_exp_diff = abs(exp_diff);
101  st_c.mantissa.hi = (st_c.mantissa.lo << 5);
102  st_c.mantissa.lo = 0;
103  uint2 cutoff_bits = (uint2)(0, 0);
104  uint2 cutoff_mask = (uint2)(0, 0);
105  if (abs_exp_diff < 32) {
106    cutoff_mask.lo = (1u << abs(exp_diff)) - 1u;
107  } else if (abs_exp_diff < 64) {
108    cutoff_mask.lo = 0xffffffff;
109    uint remaining = abs_exp_diff - 32;
110    cutoff_mask.hi = (1u << remaining) - 1u;
111  } else {
112    cutoff_mask = (uint2)(0, 0);
113  }
114  uint2 tmp = (exp_diff > 0) ? st_c.mantissa : st_mul.mantissa;
115  if (abs_exp_diff > 0) {
116    cutoff_bits = abs_exp_diff >= 64 ? tmp : (tmp & cutoff_mask);
117    if (abs_exp_diff < 32) {
118      // shift some of the hi bits into the shifted lo bits.
119      uint shift_mask = (1u << abs_exp_diff) - 1;
120      uint upper_saved_bits = tmp.hi & shift_mask;
121      upper_saved_bits = upper_saved_bits << (32 - abs_exp_diff);
122      tmp.hi >>= abs_exp_diff;
123      tmp.lo >>= abs_exp_diff;
124      tmp.lo |= upper_saved_bits;
125    } else if (abs_exp_diff < 64) {
126      tmp.lo = (tmp.hi >> (abs_exp_diff - 32));
127      tmp.hi = 0;
128    } else {
129      tmp = (uint2)(0, 0);
130    }
131  }
132  if (exp_diff > 0)
133    st_c.mantissa = tmp;
134  else
135    st_mul.mantissa = tmp;
136
137  struct fp st_fma;
138  st_fma.sign = st_mul.sign;
139  st_fma.exponent = max(st_mul.exponent, st_c.exponent);
140  st_fma.mantissa = (uint2)(0, 0);
141  if (st_c.sign == st_mul.sign) {
142    uint carry = (hadd(st_mul.mantissa.lo, st_c.mantissa.lo) >> 31) & 0x1;
143    st_fma.mantissa = st_mul.mantissa + st_c.mantissa;
144    st_fma.mantissa.hi += carry;
145  } else {
146    // cutoff bits borrow one
147    uint cutoff_borrow = ((cutoff_bits.lo != 0 || cutoff_bits.hi != 0) &&
148                          (st_mul.exponent > st_c.exponent))
149                             ? 1
150                             : 0;
151    uint borrow = 0;
152    if (st_c.mantissa.lo > st_mul.mantissa.lo) {
153      borrow = 1;
154    } else if (st_c.mantissa.lo == UINT_MAX && cutoff_borrow == 1) {
155      borrow = 1;
156    } else if ((st_c.mantissa.lo + cutoff_borrow) > st_mul.mantissa.lo) {
157      borrow = 1;
158    }
159
160    st_fma.mantissa.lo = st_mul.mantissa.lo - st_c.mantissa.lo - cutoff_borrow;
161    st_fma.mantissa.hi = st_mul.mantissa.hi - st_c.mantissa.hi - borrow;
162  }
163
164  // underflow: st_c.sign != st_mul.sign, and magnitude switches the sign
165  if (st_fma.mantissa.hi > INT_MAX) {
166    st_fma.mantissa = ~st_fma.mantissa;
167    uint carry = (hadd(st_fma.mantissa.lo, 1u) >> 31) & 0x1;
168    st_fma.mantissa.lo += 1;
169    st_fma.mantissa.hi += carry;
170
171    st_fma.sign = st_mul.sign ^ 0x80000000;
172  }
173
174  // detect overflow/underflow
175  uint leading_zeroes = clz(st_fma.mantissa.hi);
176  if (leading_zeroes == 32) {
177    leading_zeroes += clz(st_fma.mantissa.lo);
178  }
179  int overflow_bits = 3 - leading_zeroes;
180
181  // adjust exponent
182  st_fma.exponent += overflow_bits;
183
184  // handle underflow
185  if (overflow_bits < 0) {
186    uint shift = -overflow_bits;
187    if (shift < 32) {
188      uint shift_mask = (1u << shift) - 1;
189      uint saved_lo_bits = (st_fma.mantissa.lo >> (32 - shift)) & shift_mask;
190      st_fma.mantissa.lo <<= shift;
191      st_fma.mantissa.hi <<= shift;
192      st_fma.mantissa.hi |= saved_lo_bits;
193    } else if (shift < 64) {
194      st_fma.mantissa.hi = (st_fma.mantissa.lo << (64 - shift));
195      st_fma.mantissa.lo = 0;
196    } else {
197      st_fma.mantissa = (uint2)(0, 0);
198    }
199
200    overflow_bits = 0;
201  }
202
203  // rounding
204  // overflow_bits is now in the range of [0, 3] making the shift greater than
205  // 32 bits.
206  uint2 trunc_mask;
207  uint trunc_shift = C_ADJUST + overflow_bits - 32;
208  trunc_mask.hi = (1u << trunc_shift) - 1;
209  trunc_mask.lo = UINT_MAX;
210  uint2 trunc_bits = st_fma.mantissa & trunc_mask;
211  trunc_bits.lo |= (cutoff_bits.hi != 0 || cutoff_bits.lo != 0) ? 1 : 0;
212  uint2 last_bit;
213  last_bit.lo = 0;
214  last_bit.hi = st_fma.mantissa.hi & (1u << trunc_shift);
215  uint grs_shift = C_ADJUST - 3 + overflow_bits - 32;
216  uint2 grs_bits;
217  grs_bits.lo = 0;
218  grs_bits.hi = 0x4u << grs_shift;
219
220  // round to nearest even
221  if ((trunc_bits.hi > grs_bits.hi ||
222       (trunc_bits.hi == grs_bits.hi && trunc_bits.lo > grs_bits.lo)) ||
223      (trunc_bits.hi == grs_bits.hi && trunc_bits.lo == grs_bits.lo &&
224       last_bit.hi != 0)) {
225    uint shift = C_ADJUST + overflow_bits - 32;
226    st_fma.mantissa.hi += 1u << shift;
227  }
228
229        // Shift mantissa back to bit 23
230  st_fma.mantissa.lo = (st_fma.mantissa.hi >> (C_ADJUST + overflow_bits - 32));
231  st_fma.mantissa.hi = 0;
232
233  // Detect rounding overflow
234  if (st_fma.mantissa.lo > 0xffffff) {
235    ++st_fma.exponent;
236    st_fma.mantissa.lo >>= 1;
237  }
238
239  if (st_fma.mantissa.lo == 0) {
240    return 0.0f;
241  }
242
243  // Flating point range limit
244  if (st_fma.exponent > 127) {
245    return as_float(as_uint(INFINITY) | st_fma.sign);
246  }
247
248  // Flush denormals
249  if (st_fma.exponent <= -127) {
250    return as_float(st_fma.sign);
251  }
252
253  return as_float(st_fma.sign | ((st_fma.exponent + 127) << 23) |
254                  ((uint)st_fma.mantissa.lo & 0x7fffff));
255}
256_CLC_TERNARY_VECTORIZE(_CLC_DEF _CLC_OVERLOAD, float, fma, float, float, float)
257