|
22 | 22 |
|
23 | 23 | #include <clc/clc.h>
|
24 | 24 |
|
| 25 | +#include "../clcmacro.h" |
25 | 26 | #include "config.h"
|
26 | 27 | #include "math.h"
|
27 |
| -#include "../clcmacro.h" |
28 | 28 |
|
29 | 29 | struct fp {
|
30 |
| - ulong mantissa; |
31 |
| - int exponent; |
32 |
| - uint sign; |
| 30 | + ulong mantissa; |
| 31 | + int exponent; |
| 32 | + uint sign; |
33 | 33 | };
|
34 | 34 |
|
35 |
| -_CLC_DEF _CLC_OVERLOAD float __clc_sw_fma(float a, float b, float c) |
36 |
| -{ |
37 |
| - /* special cases */ |
38 |
| - if (isnan(a) || isnan(b) || isnan(c) || isinf(a) || isinf(b)) |
39 |
| - return mad(a, b, c); |
| 35 | +_CLC_DEF _CLC_OVERLOAD float __clc_sw_fma(float a, float b, float c) { |
| 36 | + /* special cases */ |
| 37 | + if (isnan(a) || isnan(b) || isnan(c) || isinf(a) || isinf(b)) |
| 38 | + return mad(a, b, c); |
40 | 39 |
|
41 |
| - /* If only c is inf, and both a,b are regular numbers, the result is c*/ |
42 |
| - if (isinf(c)) |
43 |
| - return c; |
| 40 | + /* If only c is inf, and both a,b are regular numbers, the result is c*/ |
| 41 | + if (isinf(c)) |
| 42 | + return c; |
44 | 43 |
|
45 |
| - a = __clc_flush_denormal_if_not_supported(a); |
46 |
| - b = __clc_flush_denormal_if_not_supported(b); |
47 |
| - c = __clc_flush_denormal_if_not_supported(c); |
| 44 | + a = __clc_flush_denormal_if_not_supported(a); |
| 45 | + b = __clc_flush_denormal_if_not_supported(b); |
| 46 | + c = __clc_flush_denormal_if_not_supported(c); |
48 | 47 |
|
49 |
| - if (c == 0) |
50 |
| - return a * b; |
| 48 | + if (c == 0) |
| 49 | + return a * b; |
51 | 50 |
|
52 |
| - struct fp st_a, st_b, st_c; |
| 51 | + struct fp st_a, st_b, st_c; |
53 | 52 |
|
54 |
| - st_a.exponent = a == .0f ? 0 : ((as_uint(a) & 0x7f800000) >> 23) - 127; |
55 |
| - st_b.exponent = b == .0f ? 0 : ((as_uint(b) & 0x7f800000) >> 23) - 127; |
56 |
| - st_c.exponent = c == .0f ? 0 : ((as_uint(c) & 0x7f800000) >> 23) - 127; |
| 53 | + st_a.exponent = a == .0f ? 0 : ((as_uint(a) & 0x7f800000) >> 23) - 127; |
| 54 | + st_b.exponent = b == .0f ? 0 : ((as_uint(b) & 0x7f800000) >> 23) - 127; |
| 55 | + st_c.exponent = c == .0f ? 0 : ((as_uint(c) & 0x7f800000) >> 23) - 127; |
57 | 56 |
|
58 |
| - st_a.mantissa = a == .0f ? 0 : (as_uint(a) & 0x7fffff) | 0x800000; |
59 |
| - st_b.mantissa = b == .0f ? 0 : (as_uint(b) & 0x7fffff) | 0x800000; |
60 |
| - st_c.mantissa = c == .0f ? 0 : (as_uint(c) & 0x7fffff) | 0x800000; |
| 57 | + st_a.mantissa = a == .0f ? 0 : (as_uint(a) & 0x7fffff) | 0x800000; |
| 58 | + st_b.mantissa = b == .0f ? 0 : (as_uint(b) & 0x7fffff) | 0x800000; |
| 59 | + st_c.mantissa = c == .0f ? 0 : (as_uint(c) & 0x7fffff) | 0x800000; |
61 | 60 |
|
62 |
| - st_a.sign = as_uint(a) & 0x80000000; |
63 |
| - st_b.sign = as_uint(b) & 0x80000000; |
64 |
| - st_c.sign = as_uint(c) & 0x80000000; |
| 61 | + st_a.sign = as_uint(a) & 0x80000000; |
| 62 | + st_b.sign = as_uint(b) & 0x80000000; |
| 63 | + st_c.sign = as_uint(c) & 0x80000000; |
65 | 64 |
|
66 |
| - // Multiplication. |
67 |
| - // Move the product to the highest bits to maximize precision |
68 |
| - // mantissa is 24 bits => product is 48 bits, 2bits non-fraction. |
69 |
| - // Add one bit for future addition overflow, |
70 |
| - // add another bit to detect subtraction underflow |
71 |
| - struct fp st_mul; |
72 |
| - st_mul.sign = st_a.sign ^ st_b.sign; |
73 |
| - st_mul.mantissa = (st_a.mantissa * st_b.mantissa) << 14ul; |
74 |
| - st_mul.exponent = st_mul.mantissa ? st_a.exponent + st_b.exponent : 0; |
| 65 | + // Multiplication. |
| 66 | + // Move the product to the highest bits to maximize precision |
| 67 | + // mantissa is 24 bits => product is 48 bits, 2bits non-fraction. |
| 68 | + // Add one bit for future addition overflow, |
| 69 | + // add another bit to detect subtraction underflow |
| 70 | + struct fp st_mul; |
| 71 | + st_mul.sign = st_a.sign ^ st_b.sign; |
| 72 | + st_mul.mantissa = (st_a.mantissa * st_b.mantissa) << 14ul; |
| 73 | + st_mul.exponent = st_mul.mantissa ? st_a.exponent + st_b.exponent : 0; |
75 | 74 |
|
76 |
| - // FIXME: Detecting a == 0 || b == 0 above crashed GCN isel |
77 |
| - if (st_mul.exponent == 0 && st_mul.mantissa == 0) |
78 |
| - return c; |
| 75 | + // FIXME: Detecting a == 0 || b == 0 above crashed GCN isel |
| 76 | + if (st_mul.exponent == 0 && st_mul.mantissa == 0) |
| 77 | + return c; |
79 | 78 |
|
80 | 79 | // Mantissa is 23 fractional bits, shift it the same way as product mantissa
|
81 | 80 | #define C_ADJUST 37ul
|
82 | 81 |
|
83 |
| - // both exponents are bias adjusted |
84 |
| - int exp_diff = st_mul.exponent - st_c.exponent; |
85 |
| - |
86 |
| - st_c.mantissa <<= C_ADJUST; |
87 |
| - ulong cutoff_bits = 0; |
88 |
| - ulong cutoff_mask = (1ul << abs(exp_diff)) - 1ul; |
89 |
| - if (exp_diff > 0) { |
90 |
| - cutoff_bits = exp_diff >= 64 ? st_c.mantissa : (st_c.mantissa & cutoff_mask); |
91 |
| - st_c.mantissa = exp_diff >= 64 ? 0 : (st_c.mantissa >> exp_diff); |
92 |
| - } else { |
93 |
| - cutoff_bits = -exp_diff >= 64 ? st_mul.mantissa : (st_mul.mantissa & cutoff_mask); |
94 |
| - st_mul.mantissa = -exp_diff >= 64 ? 0 : (st_mul.mantissa >> -exp_diff); |
95 |
| - } |
96 |
| - |
97 |
| - struct fp st_fma; |
98 |
| - st_fma.sign = st_mul.sign; |
99 |
| - st_fma.exponent = max(st_mul.exponent, st_c.exponent); |
100 |
| - if (st_c.sign == st_mul.sign) { |
101 |
| - st_fma.mantissa = st_mul.mantissa + st_c.mantissa; |
102 |
| - } else { |
103 |
| - // cutoff bits borrow one |
104 |
| - st_fma.mantissa = st_mul.mantissa - st_c.mantissa - (cutoff_bits && (st_mul.exponent > st_c.exponent) ? 1 : 0); |
105 |
| - } |
106 |
| - |
107 |
| - // underflow: st_c.sign != st_mul.sign, and magnitude switches the sign |
108 |
| - if (st_fma.mantissa > LONG_MAX) { |
109 |
| - st_fma.mantissa = 0 - st_fma.mantissa; |
110 |
| - st_fma.sign = st_mul.sign ^ 0x80000000; |
111 |
| - } |
112 |
| - |
113 |
| - // detect overflow/underflow |
114 |
| - int overflow_bits = 3 - clz(st_fma.mantissa); |
115 |
| - |
116 |
| - // adjust exponent |
117 |
| - st_fma.exponent += overflow_bits; |
118 |
| - |
119 |
| - // handle underflow |
120 |
| - if (overflow_bits < 0) { |
121 |
| - st_fma.mantissa <<= -overflow_bits; |
122 |
| - overflow_bits = 0; |
123 |
| - } |
124 |
| - |
125 |
| - // rounding |
126 |
| - ulong trunc_mask = (1ul << (C_ADJUST + overflow_bits)) - 1; |
127 |
| - ulong trunc_bits = (st_fma.mantissa & trunc_mask) | (cutoff_bits != 0); |
128 |
| - ulong last_bit = st_fma.mantissa & (1ul << (C_ADJUST + overflow_bits)); |
129 |
| - ulong grs_bits = (0x4ul << (C_ADJUST - 3 + overflow_bits)); |
130 |
| - |
131 |
| - // round to nearest even |
132 |
| - if ((trunc_bits > grs_bits) || |
133 |
| - (trunc_bits == grs_bits && last_bit != 0)) |
134 |
| - st_fma.mantissa += (1ul << (C_ADJUST + overflow_bits)); |
135 |
| - |
136 |
| - // Shift mantissa back to bit 23 |
137 |
| - st_fma.mantissa = (st_fma.mantissa >> (C_ADJUST + overflow_bits)); |
138 |
| - |
139 |
| - // Detect rounding overflow |
140 |
| - if (st_fma.mantissa > 0xffffff) { |
141 |
| - ++st_fma.exponent; |
142 |
| - st_fma.mantissa >>= 1; |
143 |
| - } |
144 |
| - |
145 |
| - if (st_fma.mantissa == 0) |
146 |
| - return .0f; |
147 |
| - |
148 |
| - // Flating point range limit |
149 |
| - if (st_fma.exponent > 127) |
150 |
| - return as_float(as_uint(INFINITY) | st_fma.sign); |
151 |
| - |
152 |
| - // Flush denormals |
153 |
| - if (st_fma.exponent <= -127) |
154 |
| - return as_float(st_fma.sign); |
155 |
| - |
156 |
| - return as_float(st_fma.sign | ((st_fma.exponent + 127) << 23) | ((uint)st_fma.mantissa & 0x7fffff)); |
| 82 | + // both exponents are bias adjusted |
| 83 | + int exp_diff = st_mul.exponent - st_c.exponent; |
| 84 | + |
| 85 | + st_c.mantissa <<= C_ADJUST; |
| 86 | + ulong cutoff_bits = 0; |
| 87 | + ulong cutoff_mask = (1ul << abs(exp_diff)) - 1ul; |
| 88 | + if (exp_diff > 0) { |
| 89 | + cutoff_bits = |
| 90 | + exp_diff >= 64 ? st_c.mantissa : (st_c.mantissa & cutoff_mask); |
| 91 | + st_c.mantissa = exp_diff >= 64 ? 0 : (st_c.mantissa >> exp_diff); |
| 92 | + } else { |
| 93 | + cutoff_bits = |
| 94 | + -exp_diff >= 64 ? st_mul.mantissa : (st_mul.mantissa & cutoff_mask); |
| 95 | + st_mul.mantissa = -exp_diff >= 64 ? 0 : (st_mul.mantissa >> -exp_diff); |
| 96 | + } |
| 97 | + |
| 98 | + struct fp st_fma; |
| 99 | + st_fma.sign = st_mul.sign; |
| 100 | + st_fma.exponent = max(st_mul.exponent, st_c.exponent); |
| 101 | + if (st_c.sign == st_mul.sign) { |
| 102 | + st_fma.mantissa = st_mul.mantissa + st_c.mantissa; |
| 103 | + } else { |
| 104 | + // cutoff bits borrow one |
| 105 | + st_fma.mantissa = |
| 106 | + st_mul.mantissa - st_c.mantissa - |
| 107 | + (cutoff_bits && (st_mul.exponent > st_c.exponent) ? 1 : 0); |
| 108 | + } |
| 109 | + |
| 110 | + // underflow: st_c.sign != st_mul.sign, and magnitude switches the sign |
| 111 | + if (st_fma.mantissa > LONG_MAX) { |
| 112 | + st_fma.mantissa = 0 - st_fma.mantissa; |
| 113 | + st_fma.sign = st_mul.sign ^ 0x80000000; |
| 114 | + } |
| 115 | + |
| 116 | + // detect overflow/underflow |
| 117 | + int overflow_bits = 3 - clz(st_fma.mantissa); |
| 118 | + |
| 119 | + // adjust exponent |
| 120 | + st_fma.exponent += overflow_bits; |
| 121 | + |
| 122 | + // handle underflow |
| 123 | + if (overflow_bits < 0) { |
| 124 | + st_fma.mantissa <<= -overflow_bits; |
| 125 | + overflow_bits = 0; |
| 126 | + } |
| 127 | + |
| 128 | + // rounding |
| 129 | + ulong trunc_mask = (1ul << (C_ADJUST + overflow_bits)) - 1; |
| 130 | + ulong trunc_bits = (st_fma.mantissa & trunc_mask) | (cutoff_bits != 0); |
| 131 | + ulong last_bit = st_fma.mantissa & (1ul << (C_ADJUST + overflow_bits)); |
| 132 | + ulong grs_bits = (0x4ul << (C_ADJUST - 3 + overflow_bits)); |
| 133 | + |
| 134 | + // round to nearest even |
| 135 | + if ((trunc_bits > grs_bits) || (trunc_bits == grs_bits && last_bit != 0)) |
| 136 | + st_fma.mantissa += (1ul << (C_ADJUST + overflow_bits)); |
| 137 | + |
| 138 | + // Shift mantissa back to bit 23 |
| 139 | + st_fma.mantissa = (st_fma.mantissa >> (C_ADJUST + overflow_bits)); |
| 140 | + |
| 141 | + // Detect rounding overflow |
| 142 | + if (st_fma.mantissa > 0xffffff) { |
| 143 | + ++st_fma.exponent; |
| 144 | + st_fma.mantissa >>= 1; |
| 145 | + } |
| 146 | + |
| 147 | + if (st_fma.mantissa == 0) |
| 148 | + return .0f; |
| 149 | + |
| 150 | + // Flating point range limit |
| 151 | + if (st_fma.exponent > 127) |
| 152 | + return as_float(as_uint(INFINITY) | st_fma.sign); |
| 153 | + |
| 154 | + // Flush denormals |
| 155 | + if (st_fma.exponent <= -127) |
| 156 | + return as_float(st_fma.sign); |
| 157 | + |
| 158 | + return as_float(st_fma.sign | ((st_fma.exponent + 127) << 23) | |
| 159 | + ((uint)st_fma.mantissa & 0x7fffff)); |
157 | 160 | }
|
158 |
| -_CLC_TERNARY_VECTORIZE(_CLC_DEF _CLC_OVERLOAD, float, __clc_sw_fma, float, float, float) |
| 161 | +_CLC_TERNARY_VECTORIZE(_CLC_DEF _CLC_OVERLOAD, float, __clc_sw_fma, float, |
| 162 | + float, float) |
0 commit comments