xref: /openbmc/qemu/fpu/softfloat-parts.c.inc (revision da10a907) 
1/*
2 * QEMU float support
3 *
4 * The code in this source file is derived from release 2a of the SoftFloat
5 * IEC/IEEE Floating-point Arithmetic Package. Those parts of the code (and
6 * some later contributions) are provided under that license, as detailed below.
7 * It has subsequently been modified by contributors to the QEMU Project,
8 * so some portions are provided under:
9 *  the SoftFloat-2a license
10 *  the BSD license
11 *  GPL-v2-or-later
12 *
13 * Any future contributions to this file after December 1st 2014 will be
14 * taken to be licensed under the Softfloat-2a license unless specifically
15 * indicated otherwise.
16 */
17
18static void partsN(return_nan)(FloatPartsN *a, float_status *s)
19{
20    switch (a->cls) {
21    case float_class_snan:
22        float_raise(float_flag_invalid, s);
23        if (s->default_nan_mode) {
24            parts_default_nan(a, s);
25        } else {
26            parts_silence_nan(a, s);
27        }
28        break;
29    case float_class_qnan:
30        if (s->default_nan_mode) {
31            parts_default_nan(a, s);
32        }
33        break;
34    default:
35        g_assert_not_reached();
36    }
37}
38
39static FloatPartsN *partsN(pick_nan)(FloatPartsN *a, FloatPartsN *b,
40                                     float_status *s)
41{
42    if (is_snan(a->cls) || is_snan(b->cls)) {
43        float_raise(float_flag_invalid, s);
44    }
45
46    if (s->default_nan_mode) {
47        parts_default_nan(a, s);
48    } else {
49        int cmp = frac_cmp(a, b);
50        if (cmp == 0) {
51            cmp = a->sign < b->sign;
52        }
53
54        if (pickNaN(a->cls, b->cls, cmp > 0, s)) {
55            a = b;
56        }
57        if (is_snan(a->cls)) {
58            parts_silence_nan(a, s);
59        }
60    }
61    return a;
62}
63
64static FloatPartsN *partsN(pick_nan_muladd)(FloatPartsN *a, FloatPartsN *b,
65                                            FloatPartsN *c, float_status *s,
66                                            int ab_mask, int abc_mask)
67{
68    int which;
69
70    if (unlikely(abc_mask & float_cmask_snan)) {
71        float_raise(float_flag_invalid, s);
72    }
73
74    which = pickNaNMulAdd(a->cls, b->cls, c->cls,
75                          ab_mask == float_cmask_infzero, s);
76
77    if (s->default_nan_mode || which == 3) {
78        /*
79         * Note that this check is after pickNaNMulAdd so that function
80         * has an opportunity to set the Invalid flag for infzero.
81         */
82        parts_default_nan(a, s);
83        return a;
84    }
85
86    switch (which) {
87    case 0:
88        break;
89    case 1:
90        a = b;
91        break;
92    case 2:
93        a = c;
94        break;
95    default:
96        g_assert_not_reached();
97    }
98    if (is_snan(a->cls)) {
99        parts_silence_nan(a, s);
100    }
101    return a;
102}
103
104/*
105 * Canonicalize the FloatParts structure.  Determine the class,
106 * unbias the exponent, and normalize the fraction.
107 */
108static void partsN(canonicalize)(FloatPartsN *p, float_status *status,
109                                 const FloatFmt *fmt)
110{
111    if (unlikely(p->exp == 0)) {
112        if (likely(frac_eqz(p))) {
113            p->cls = float_class_zero;
114        } else if (status->flush_inputs_to_zero) {
115            float_raise(float_flag_input_denormal, status);
116            p->cls = float_class_zero;
117            frac_clear(p);
118        } else {
119            int shift = frac_normalize(p);
120            p->cls = float_class_normal;
121            p->exp = fmt->frac_shift - fmt->exp_bias - shift + 1;
122        }
123    } else if (likely(p->exp < fmt->exp_max) || fmt->arm_althp) {
124        p->cls = float_class_normal;
125        p->exp -= fmt->exp_bias;
126        frac_shl(p, fmt->frac_shift);
127        p->frac_hi |= DECOMPOSED_IMPLICIT_BIT;
128    } else if (likely(frac_eqz(p))) {
129        p->cls = float_class_inf;
130    } else {
131        frac_shl(p, fmt->frac_shift);
132        p->cls = (parts_is_snan_frac(p->frac_hi, status)
133                  ? float_class_snan : float_class_qnan);
134    }
135}
136
137/*
138 * Round and uncanonicalize a floating-point number by parts. There
139 * are FRAC_SHIFT bits that may require rounding at the bottom of the
140 * fraction; these bits will be removed. The exponent will be biased
141 * by EXP_BIAS and must be bounded by [EXP_MAX-1, 0].
142 */
143static void partsN(uncanon)(FloatPartsN *p, float_status *s,
144                            const FloatFmt *fmt)
145{
146    const int exp_max = fmt->exp_max;
147    const int frac_shift = fmt->frac_shift;
148    const uint64_t frac_lsb = fmt->frac_lsb;
149    const uint64_t frac_lsbm1 = fmt->frac_lsbm1;
150    const uint64_t round_mask = fmt->round_mask;
151    const uint64_t roundeven_mask = fmt->roundeven_mask;
152    uint64_t inc;
153    bool overflow_norm;
154    int exp, flags = 0;
155
156    if (unlikely(p->cls != float_class_normal)) {
157        switch (p->cls) {
158        case float_class_zero:
159            p->exp = 0;
160            frac_clear(p);
161            return;
162        case float_class_inf:
163            g_assert(!fmt->arm_althp);
164            p->exp = fmt->exp_max;
165            frac_clear(p);
166            return;
167        case float_class_qnan:
168        case float_class_snan:
169            g_assert(!fmt->arm_althp);
170            p->exp = fmt->exp_max;
171            frac_shr(p, fmt->frac_shift);
172            return;
173        default:
174            break;
175        }
176        g_assert_not_reached();
177    }
178
179    switch (s->float_rounding_mode) {
180    case float_round_nearest_even:
181        overflow_norm = false;
182        inc = ((p->frac_lo & roundeven_mask) != frac_lsbm1 ? frac_lsbm1 : 0);
183        break;
184    case float_round_ties_away:
185        overflow_norm = false;
186        inc = frac_lsbm1;
187        break;
188    case float_round_to_zero:
189        overflow_norm = true;
190        inc = 0;
191        break;
192    case float_round_up:
193        inc = p->sign ? 0 : round_mask;
194        overflow_norm = p->sign;
195        break;
196    case float_round_down:
197        inc = p->sign ? round_mask : 0;
198        overflow_norm = !p->sign;
199        break;
200    case float_round_to_odd:
201        overflow_norm = true;
202        inc = p->frac_lo & frac_lsb ? 0 : round_mask;
203        break;
204    default:
205        g_assert_not_reached();
206    }
207
208    exp = p->exp + fmt->exp_bias;
209    if (likely(exp > 0)) {
210        if (p->frac_lo & round_mask) {
211            flags |= float_flag_inexact;
212            if (frac_addi(p, p, inc)) {
213                frac_shr(p, 1);
214                p->frac_hi |= DECOMPOSED_IMPLICIT_BIT;
215                exp++;
216            }
217        }
218        frac_shr(p, frac_shift);
219
220        if (fmt->arm_althp) {
221            /* ARM Alt HP eschews Inf and NaN for a wider exponent.  */
222            if (unlikely(exp > exp_max)) {
223                /* Overflow.  Return the maximum normal.  */
224                flags = float_flag_invalid;
225                exp = exp_max;
226                frac_allones(p);
227            }
228        } else if (unlikely(exp >= exp_max)) {
229            flags |= float_flag_overflow | float_flag_inexact;
230            if (overflow_norm) {
231                exp = exp_max - 1;
232                frac_allones(p);
233            } else {
234                p->cls = float_class_inf;
235                exp = exp_max;
236                frac_clear(p);
237            }
238        }
239    } else if (s->flush_to_zero) {
240        flags |= float_flag_output_denormal;
241        p->cls = float_class_zero;
242        exp = 0;
243        frac_clear(p);
244    } else {
245        bool is_tiny = s->tininess_before_rounding || exp < 0;
246
247        if (!is_tiny) {
248            FloatPartsN discard;
249            is_tiny = !frac_addi(&discard, p, inc);
250        }
251
252        frac_shrjam(p, 1 - exp);
253
254        if (p->frac_lo & round_mask) {
255            /* Need to recompute round-to-even/round-to-odd. */
256            switch (s->float_rounding_mode) {
257            case float_round_nearest_even:
258                inc = ((p->frac_lo & roundeven_mask) != frac_lsbm1
259                       ? frac_lsbm1 : 0);
260                break;
261            case float_round_to_odd:
262                inc = p->frac_lo & frac_lsb ? 0 : round_mask;
263                break;
264            default:
265                break;
266            }
267            flags |= float_flag_inexact;
268            frac_addi(p, p, inc);
269        }
270
271        exp = (p->frac_hi & DECOMPOSED_IMPLICIT_BIT) != 0;
272        frac_shr(p, frac_shift);
273
274        if (is_tiny && (flags & float_flag_inexact)) {
275            flags |= float_flag_underflow;
276        }
277        if (exp == 0 && frac_eqz(p)) {
278            p->cls = float_class_zero;
279        }
280    }
281    p->exp = exp;
282    float_raise(flags, s);
283}
284
285/*
286 * Returns the result of adding or subtracting the values of the
287 * floating-point values `a' and `b'. The operation is performed
288 * according to the IEC/IEEE Standard for Binary Floating-Point
289 * Arithmetic.
290 */
291static FloatPartsN *partsN(addsub)(FloatPartsN *a, FloatPartsN *b,
292                                   float_status *s, bool subtract)
293{
294    bool b_sign = b->sign ^ subtract;
295    int ab_mask = float_cmask(a->cls) | float_cmask(b->cls);
296
297    if (a->sign != b_sign) {
298        /* Subtraction */
299        if (likely(ab_mask == float_cmask_normal)) {
300            if (parts_sub_normal(a, b)) {
301                return a;
302            }
303            /* Subtract was exact, fall through to set sign. */
304            ab_mask = float_cmask_zero;
305        }
306
307        if (ab_mask == float_cmask_zero) {
308            a->sign = s->float_rounding_mode == float_round_down;
309            return a;
310        }
311
312        if (unlikely(ab_mask & float_cmask_anynan)) {
313            goto p_nan;
314        }
315
316        if (ab_mask & float_cmask_inf) {
317            if (a->cls != float_class_inf) {
318                /* N - Inf */
319                goto return_b;
320            }
321            if (b->cls != float_class_inf) {
322                /* Inf - N */
323                return a;
324            }
325            /* Inf - Inf */
326            float_raise(float_flag_invalid, s);
327            parts_default_nan(a, s);
328            return a;
329        }
330    } else {
331        /* Addition */
332        if (likely(ab_mask == float_cmask_normal)) {
333            parts_add_normal(a, b);
334            return a;
335        }
336
337        if (ab_mask == float_cmask_zero) {
338            return a;
339        }
340
341        if (unlikely(ab_mask & float_cmask_anynan)) {
342            goto p_nan;
343        }
344
345        if (ab_mask & float_cmask_inf) {
346            a->cls = float_class_inf;
347            return a;
348        }
349    }
350
351    if (b->cls == float_class_zero) {
352        g_assert(a->cls == float_class_normal);
353        return a;
354    }
355
356    g_assert(a->cls == float_class_zero);
357    g_assert(b->cls == float_class_normal);
358 return_b:
359    b->sign = b_sign;
360    return b;
361
362 p_nan:
363    return parts_pick_nan(a, b, s);
364}
365