fplib.cc revision 10272:336c7d36ac00 
1/*
2* Copyright (c) 2012-2013 ARM Limited
3* All rights reserved
4*
5* The license below extends only to copyright in the software and shall
6* not be construed as granting a license to any other intellectual
7* property including but not limited to intellectual property relating
8* to a hardware implementation of the functionality of the software
9* licensed hereunder.  You may use the software subject to the license
10* terms below provided that you ensure that this notice is replicated
11* unmodified and in its entirety in all distributions of the software,
12* modified or unmodified, in source code or in binary form.
13*
14* Redistribution and use in source and binary forms, with or without
15* modification, are permitted provided that the following conditions are
16* met: redistributions of source code must retain the above copyright
17* notice, this list of conditions and the following disclaimer;
18* redistributions in binary form must reproduce the above copyright
19* notice, this list of conditions and the following disclaimer in the
20* documentation and/or other materials provided with the distribution;
21* neither the name of the copyright holders nor the names of its
22* contributors may be used to endorse or promote products derived from
23* this software without specific prior written permission.
24*
25* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
26* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
27* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
28* A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
29* OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
30* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
31* LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
32* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
33* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
34* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
35* OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36*
37* Authors: Edmund Grimley Evans
38*          Thomas Grocutt
39*/
40
41#include <stdint.h>
42
43#include <cassert>
44
45#include "fplib.hh"
46
47namespace ArmISA
48{
49
50#define FPLIB_RN 0
51#define FPLIB_RP 1
52#define FPLIB_RM 2
53#define FPLIB_RZ 3
54#define FPLIB_FZ 4
55#define FPLIB_DN 8
56#define FPLIB_AHP 16
57
58#define FPLIB_IDC 128 // Input Denormal
59#define FPLIB_IXC 16  // Inexact
60#define FPLIB_UFC 8   // Underflow
61#define FPLIB_OFC 4   // Overflow
62#define FPLIB_DZC 2   // Division by Zero
63#define FPLIB_IOC 1   // Invalid Operation
64
65static inline uint16_t
66lsl16(uint16_t x, uint32_t shift)
67{
68    return shift < 16 ? x << shift : 0;
69}
70
71static inline uint16_t
72lsr16(uint16_t x, uint32_t shift)
73{
74    return shift < 16 ? x >> shift : 0;
75}
76
77static inline uint32_t
78lsl32(uint32_t x, uint32_t shift)
79{
80    return shift < 32 ? x << shift : 0;
81}
82
83static inline uint32_t
84lsr32(uint32_t x, uint32_t shift)
85{
86    return shift < 32 ? x >> shift : 0;
87}
88
89static inline uint64_t
90lsl64(uint64_t x, uint32_t shift)
91{
92    return shift < 64 ? x << shift : 0;
93}
94
95static inline uint64_t
96lsr64(uint64_t x, uint32_t shift)
97{
98    return shift < 64 ? x >> shift : 0;
99}
100
101static inline void
102lsl128(uint64_t *r0, uint64_t *r1, uint64_t x0, uint64_t x1, uint32_t shift)
103{
104    if (shift < 64) {
105        *r1 = x1 << shift | x0 >> (64 - shift);
106        *r0 = x0 << shift;
107    } else if (shift < 128) {
108        *r1 = x0 << (shift - 64);
109        *r0 = 0;
110    } else {
111        *r1 = 0;
112        *r0 = 0;
113    }
114}
115
116static inline void
117lsr128(uint64_t *r0, uint64_t *r1, uint64_t x0, uint64_t x1, uint32_t shift)
118{
119    if (shift < 64) {
120        *r0 = x0 >> shift | x1 << (64 - shift);
121        *r1 = x1 >> shift;
122    } else if (shift < 128) {
123        *r0 = x1 >> (shift - 64);
124        *r1 = 0;
125    } else {
126        *r0 = 0;
127        *r1 = 0;
128    }
129}
130
131static inline void
132mul62x62(uint64_t *x0, uint64_t *x1, uint64_t a, uint64_t b)
133{
134    uint32_t mask = ((uint32_t)1 << 31) - 1;
135    uint64_t a0 = a & mask;
136    uint64_t a1 = a >> 31 & mask;
137    uint64_t b0 = b & mask;
138    uint64_t b1 = b >> 31 & mask;
139    uint64_t p0 = a0 * b0;
140    uint64_t p2 = a1 * b1;
141    uint64_t p1 = (a0 + a1) * (b0 + b1) - p0 - p2;
142    uint64_t s0 = p0;
143    uint64_t s1 = (s0 >> 31) + p1;
144    uint64_t s2 = (s1 >> 31) + p2;
145    *x0 = (s0 & mask) | (s1 & mask) << 31 | s2 << 62;
146    *x1 = s2 >> 2;
147}
148
149static inline
150void mul64x32(uint64_t *x0, uint64_t *x1, uint64_t a, uint32_t b)
151{
152    uint64_t t0 = (uint64_t)(uint32_t)a * b;
153    uint64_t t1 = (t0 >> 32) + (a >> 32) * b;
154    *x0 = t1 << 32 | (uint32_t)t0;
155    *x1 = t1 >> 32;
156}
157
158static inline void
159add128(uint64_t *x0, uint64_t *x1, uint64_t a0, uint64_t a1, uint64_t b0,
160       uint64_t b1)
161{
162    *x0 = a0 + b0;
163    *x1 = a1 + b1 + (*x0 < a0);
164}
165
166static inline void
167sub128(uint64_t *x0, uint64_t *x1, uint64_t a0, uint64_t a1, uint64_t b0,
168       uint64_t b1)
169{
170    *x0 = a0 - b0;
171    *x1 = a1 - b1 - (*x0 > a0);
172}
173
174static inline int
175cmp128(uint64_t a0, uint64_t a1, uint64_t b0, uint64_t b1)
176{
177    return (a1 < b1 ? -1 : a1 > b1 ? 1 : a0 < b0 ? -1 : a0 > b0 ? 1 : 0);
178}
179
180static inline uint16_t
181fp16_normalise(uint16_t mnt, int *exp)
182{
183    int shift;
184
185    if (!mnt) {
186        return 0;
187    }
188
189    for (shift = 8; shift; shift >>= 1) {
190        if (!(mnt >> (16 - shift))) {
191            mnt <<= shift;
192            *exp -= shift;
193        }
194    }
195    return mnt;
196}
197
198static inline uint32_t
199fp32_normalise(uint32_t mnt, int *exp)
200{
201    int shift;
202
203    if (!mnt) {
204        return 0;
205    }
206
207    for (shift = 16; shift; shift >>= 1) {
208        if (!(mnt >> (32 - shift))) {
209            mnt <<= shift;
210            *exp -= shift;
211        }
212    }
213    return mnt;
214}
215
216static inline uint64_t
217fp64_normalise(uint64_t mnt, int *exp)
218{
219    int shift;
220
221    if (!mnt) {
222        return 0;
223    }
224
225    for (shift = 32; shift; shift >>= 1) {
226        if (!(mnt >> (64 - shift))) {
227            mnt <<= shift;
228            *exp -= shift;
229        }
230    }
231    return mnt;
232}
233
234static inline void
235fp128_normalise(uint64_t *mnt0, uint64_t *mnt1, int *exp)
236{
237    uint64_t x0 = *mnt0;
238    uint64_t x1 = *mnt1;
239    int shift;
240
241    if (!x0 && !x1) {
242        return;
243    }
244
245    if (!x1) {
246        x1 = x0;
247        x0 = 0;
248        *exp -= 64;
249    }
250
251    for (shift = 32; shift; shift >>= 1) {
252        if (!(x1 >> (64 - shift))) {
253            x1 = x1 << shift | x0 >> (64 - shift);
254            x0 <<= shift;
255            *exp -= shift;
256        }
257    }
258
259    *mnt0 = x0;
260    *mnt1 = x1;
261}
262
263static inline uint16_t
264fp16_pack(uint16_t sgn, uint16_t exp, uint16_t mnt)
265{
266    return sgn << 15 | exp << 10 | (mnt & (((uint16_t)1 << 10) - 1));
267}
268
269static inline uint32_t
270fp32_pack(uint32_t sgn, uint32_t exp, uint32_t mnt)
271{
272    return sgn << 31 | exp << 23 | (mnt & (((uint32_t)1 << 23) - 1));
273}
274
275static inline uint64_t
276fp64_pack(uint64_t sgn, uint64_t exp, uint64_t mnt)
277{
278    return (uint64_t)sgn << 63 | exp << 52 | (mnt & (((uint64_t)1 << 52) - 1));
279}
280
281static inline uint16_t
282fp16_zero(int sgn)
283{
284    return fp16_pack(sgn, 0, 0);
285}
286
287static inline uint32_t
288fp32_zero(int sgn)
289{
290    return fp32_pack(sgn, 0, 0);
291}
292
293static inline uint64_t
294fp64_zero(int sgn)
295{
296    return fp64_pack(sgn, 0, 0);
297}
298
299static inline uint16_t
300fp16_max_normal(int sgn)
301{
302    return fp16_pack(sgn, 30, -1);
303}
304
305static inline uint32_t
306fp32_max_normal(int sgn)
307{
308    return fp32_pack(sgn, 254, -1);
309}
310
311static inline uint64_t
312fp64_max_normal(int sgn)
313{
314    return fp64_pack(sgn, 2046, -1);
315}
316
317static inline uint16_t
318fp16_infinity(int sgn)
319{
320    return fp16_pack(sgn, 31, 0);
321}
322
323static inline uint32_t
324fp32_infinity(int sgn)
325{
326    return fp32_pack(sgn, 255, 0);
327}
328
329static inline uint64_t
330fp64_infinity(int sgn)
331{
332    return fp64_pack(sgn, 2047, 0);
333}
334
335static inline uint16_t
336fp16_defaultNaN()
337{
338    return fp16_pack(0, 31, (uint16_t)1 << 9);
339}
340
341static inline uint32_t
342fp32_defaultNaN()
343{
344    return fp32_pack(0, 255, (uint32_t)1 << 22);
345}
346
347static inline uint64_t
348fp64_defaultNaN()
349{
350    return fp64_pack(0, 2047, (uint64_t)1 << 51);
351}
352
353static inline void
354fp16_unpack(int *sgn, int *exp, uint16_t *mnt, uint16_t x, int mode,
355            int *flags)
356{
357    *sgn = x >> 15;
358    *exp = x >> 10 & 31;
359    *mnt = x & (((uint16_t)1 << 10) - 1);
360
361    // Handle subnormals:
362    if (*exp) {
363        *mnt |= (uint16_t)1 << 10;
364    } else {
365        ++*exp;
366        // There is no flush to zero in this case!
367    }
368}
369
370static inline void
371fp32_unpack(int *sgn, int *exp, uint32_t *mnt, uint32_t x, int mode,
372            int *flags)
373{
374    *sgn = x >> 31;
375    *exp = x >> 23 & 255;
376    *mnt = x & (((uint32_t)1 << 23) - 1);
377
378    // Handle subnormals:
379    if (*exp) {
380        *mnt |= (uint32_t)1 << 23;
381    } else {
382        ++*exp;
383        if ((mode & FPLIB_FZ) && *mnt) {
384            *flags |= FPLIB_IDC;
385            *mnt = 0;
386        }
387    }
388}
389
390static inline void
391fp64_unpack(int *sgn, int *exp, uint64_t *mnt, uint64_t x, int mode,
392            int *flags)
393{
394    *sgn = x >> 63;
395    *exp = x >> 52 & 2047;
396    *mnt = x & (((uint64_t)1 << 52) - 1);
397
398    // Handle subnormals:
399    if (*exp) {
400        *mnt |= (uint64_t)1 << 52;
401    } else {
402        ++*exp;
403        if ((mode & FPLIB_FZ) && *mnt) {
404            *flags |= FPLIB_IDC;
405            *mnt = 0;
406        }
407    }
408}
409
410static inline uint32_t
411fp32_process_NaN(uint32_t a, int mode, int *flags)
412{
413    if (!(a >> 22 & 1)) {
414        *flags |= FPLIB_IOC;
415        a |= (uint32_t)1 << 22;
416    }
417    return mode & FPLIB_DN ? fp32_defaultNaN() : a;
418}
419
420static inline uint64_t
421fp64_process_NaN(uint64_t a, int mode, int *flags)
422{
423    if (!(a >> 51 & 1)) {
424        *flags |= FPLIB_IOC;
425        a |= (uint64_t)1 << 51;
426    }
427    return mode & FPLIB_DN ? fp64_defaultNaN() : a;
428}
429
430static uint32_t
431fp32_process_NaNs(uint32_t a, uint32_t b, int mode, int *flags)
432{
433    int a_exp = a >> 23 & 255;
434    uint32_t a_mnt = a & (((uint32_t)1 << 23) - 1);
435    int b_exp = b >> 23 & 255;
436    uint32_t b_mnt = b & (((uint32_t)1 << 23) - 1);
437
438    // Handle signalling NaNs:
439    if (a_exp == 255 && a_mnt && !(a_mnt >> 22 & 1))
440        return fp32_process_NaN(a, mode, flags);
441    if (b_exp == 255 && b_mnt && !(b_mnt >> 22 & 1))
442        return fp32_process_NaN(b, mode, flags);
443
444    // Handle quiet NaNs:
445    if (a_exp == 255 && a_mnt)
446        return fp32_process_NaN(a, mode, flags);
447    if (b_exp == 255 && b_mnt)
448        return fp32_process_NaN(b, mode, flags);
449
450    return 0;
451}
452
453static uint64_t
454fp64_process_NaNs(uint64_t a, uint64_t b, int mode, int *flags)
455{
456    int a_exp = a >> 52 & 2047;
457    uint64_t a_mnt = a & (((uint64_t)1 << 52) - 1);
458    int b_exp = b >> 52 & 2047;
459    uint64_t b_mnt = b & (((uint64_t)1 << 52) - 1);
460
461    // Handle signalling NaNs:
462    if (a_exp == 2047 && a_mnt && !(a_mnt >> 51 & 1))
463        return fp64_process_NaN(a, mode, flags);
464    if (b_exp == 2047 && b_mnt && !(b_mnt >> 51 & 1))
465        return fp64_process_NaN(b, mode, flags);
466
467    // Handle quiet NaNs:
468    if (a_exp == 2047 && a_mnt)
469        return fp64_process_NaN(a, mode, flags);
470    if (b_exp == 2047 && b_mnt)
471        return fp64_process_NaN(b, mode, flags);
472
473    return 0;
474}
475
476static uint32_t
477fp32_process_NaNs3(uint32_t a, uint32_t b, uint32_t c, int mode, int *flags)
478{
479    int a_exp = a >> 23 & 255;
480    uint32_t a_mnt = a & (((uint32_t)1 << 23) - 1);
481    int b_exp = b >> 23 & 255;
482    uint32_t b_mnt = b & (((uint32_t)1 << 23) - 1);
483    int c_exp = c >> 23 & 255;
484    uint32_t c_mnt = c & (((uint32_t)1 << 23) - 1);
485
486    // Handle signalling NaNs:
487    if (a_exp == 255 && a_mnt && !(a_mnt >> 22 & 1))
488        return fp32_process_NaN(a, mode, flags);
489    if (b_exp == 255 && b_mnt && !(b_mnt >> 22 & 1))
490        return fp32_process_NaN(b, mode, flags);
491    if (c_exp == 255 && c_mnt && !(c_mnt >> 22 & 1))
492        return fp32_process_NaN(c, mode, flags);
493
494    // Handle quiet NaNs:
495    if (a_exp == 255 && a_mnt)
496        return fp32_process_NaN(a, mode, flags);
497    if (b_exp == 255 && b_mnt)
498        return fp32_process_NaN(b, mode, flags);
499    if (c_exp == 255 && c_mnt)
500        return fp32_process_NaN(c, mode, flags);
501
502    return 0;
503}
504
505static uint64_t
506fp64_process_NaNs3(uint64_t a, uint64_t b, uint64_t c, int mode, int *flags)
507{
508    int a_exp = a >> 52 & 2047;
509    uint64_t a_mnt = a & (((uint64_t)1 << 52) - 1);
510    int b_exp = b >> 52 & 2047;
511    uint64_t b_mnt = b & (((uint64_t)1 << 52) - 1);
512    int c_exp = c >> 52 & 2047;
513    uint64_t c_mnt = c & (((uint64_t)1 << 52) - 1);
514
515    // Handle signalling NaNs:
516    if (a_exp == 2047 && a_mnt && !(a_mnt >> 51 & 1))
517        return fp64_process_NaN(a, mode, flags);
518    if (b_exp == 2047 && b_mnt && !(b_mnt >> 51 & 1))
519        return fp64_process_NaN(b, mode, flags);
520    if (c_exp == 2047 && c_mnt && !(c_mnt >> 51 & 1))
521        return fp64_process_NaN(c, mode, flags);
522
523    // Handle quiet NaNs:
524    if (a_exp == 2047 && a_mnt)
525        return fp64_process_NaN(a, mode, flags);
526    if (b_exp == 2047 && b_mnt)
527        return fp64_process_NaN(b, mode, flags);
528    if (c_exp == 2047 && c_mnt)
529        return fp64_process_NaN(c, mode, flags);
530
531    return 0;
532}
533
534static uint16_t
535fp16_round_(int sgn, int exp, uint16_t mnt, int rm, int mode, int *flags)
536{
537    int biased_exp; // non-negative exponent value for result
538    uint16_t int_mant; // mantissa for result, less than (1 << 11)
539    int error; // 0, 1, 2 or 3, where 2 means int_mant is wrong by exactly 0.5
540
541    assert(rm != FPRounding_TIEAWAY);
542
543    // There is no flush to zero in this case!
544
545    // The bottom 5 bits of mnt are orred together:
546    mnt = (uint16_t)1 << 12 | mnt >> 4 | ((mnt & 31) != 0);
547
548    if (exp > 0) {
549        biased_exp = exp;
550        int_mant = mnt >> 2;
551        error = mnt & 3;
552    } else {
553        biased_exp = 0;
554        int_mant = lsr16(mnt, 3 - exp);
555        error = (lsr16(mnt, 1 - exp) & 3) | !!(mnt & (lsl16(1, 1 - exp) - 1));
556    }
557
558    if (!biased_exp && error) { // xx should also check fpscr_val<11>
559        *flags |= FPLIB_UFC;
560    }
561
562    // Round up:
563    if ((rm == FPLIB_RN && (error == 3 ||
564                            (error == 2 && (int_mant & 1)))) ||
565        (((rm == FPLIB_RP && !sgn) || (rm == FPLIB_RM && sgn)) && error)) {
566        ++int_mant;
567        if (int_mant == (uint32_t)1 << 10) {
568            // Rounded up from denormalized to normalized
569            biased_exp = 1;
570        }
571        if (int_mant == (uint32_t)1 << 11) {
572            // Rounded up to next exponent
573            ++biased_exp;
574            int_mant >>= 1;
575        }
576    }
577
578    // Handle rounding to odd aka Von Neumann rounding:
579    if (error && rm == FPRounding_ODD)
580        int_mant |= 1;
581
582    // Handle overflow:
583    if (!(mode & FPLIB_AHP)) {
584        if (biased_exp >= 31) {
585            *flags |= FPLIB_OFC | FPLIB_IXC;
586            if (rm == FPLIB_RN || (rm == FPLIB_RP && !sgn) ||
587                (rm == FPLIB_RM && sgn)) {
588                return fp16_infinity(sgn);
589            } else {
590                return fp16_max_normal(sgn);
591            }
592        }
593    } else {
594        if (biased_exp >= 32) {
595            *flags |= FPLIB_IOC;
596            return fp16_pack(sgn, 31, -1);
597        }
598    }
599
600    if (error) {
601        *flags |= FPLIB_IXC;
602    }
603
604    return fp16_pack(sgn, biased_exp, int_mant);
605}
606
607static uint32_t
608fp32_round_(int sgn, int exp, uint32_t mnt, int rm, int mode, int *flags)
609{
610    int biased_exp; // non-negative exponent value for result
611    uint32_t int_mant; // mantissa for result, less than (1 << 24)
612    int error; // 0, 1, 2 or 3, where 2 means int_mant is wrong by exactly 0.5
613
614    assert(rm != FPRounding_TIEAWAY);
615
616    // Flush to zero:
617    if ((mode & FPLIB_FZ) && exp < 1) {
618        *flags |= FPLIB_UFC;
619        return fp32_zero(sgn);
620    }
621
622    // The bottom 8 bits of mnt are orred together:
623    mnt = (uint32_t)1 << 25 | mnt >> 7 | ((mnt & 255) != 0);
624
625    if (exp > 0) {
626        biased_exp = exp;
627        int_mant = mnt >> 2;
628        error = mnt & 3;
629    } else {
630        biased_exp = 0;
631        int_mant = lsr32(mnt, 3 - exp);
632        error = (lsr32(mnt, 1 - exp) & 3) | !!(mnt & (lsl32(1, 1 - exp) - 1));
633    }
634
635    if (!biased_exp && error) { // xx should also check fpscr_val<11>
636        *flags |= FPLIB_UFC;
637    }
638
639    // Round up:
640    if ((rm == FPLIB_RN && (error == 3 ||
641                            (error == 2 && (int_mant & 1)))) ||
642        (((rm == FPLIB_RP && !sgn) || (rm == FPLIB_RM && sgn)) && error)) {
643        ++int_mant;
644        if (int_mant == (uint32_t)1 << 23) {
645            // Rounded up from denormalized to normalized
646            biased_exp = 1;
647        }
648        if (int_mant == (uint32_t)1 << 24) {
649            // Rounded up to next exponent
650            ++biased_exp;
651            int_mant >>= 1;
652        }
653    }
654
655    // Handle rounding to odd aka Von Neumann rounding:
656    if (error && rm == FPRounding_ODD)
657        int_mant |= 1;
658
659    // Handle overflow:
660    if (biased_exp >= 255) {
661        *flags |= FPLIB_OFC | FPLIB_IXC;
662        if (rm == FPLIB_RN || (rm == FPLIB_RP && !sgn) ||
663            (rm == FPLIB_RM && sgn)) {
664            return fp32_infinity(sgn);
665        } else {
666            return fp32_max_normal(sgn);
667        }
668    }
669
670    if (error) {
671        *flags |= FPLIB_IXC;
672    }
673
674    return fp32_pack(sgn, biased_exp, int_mant);
675}
676
677static uint32_t
678fp32_round(int sgn, int exp, uint32_t mnt, int mode, int *flags)
679{
680    return fp32_round_(sgn, exp, mnt, mode & 3, mode, flags);
681}
682
683static uint64_t
684fp64_round_(int sgn, int exp, uint64_t mnt, int rm, int mode, int *flags)
685{
686    int biased_exp; // non-negative exponent value for result
687    uint64_t int_mant; // mantissa for result, less than (1 << 52)
688    int error; // 0, 1, 2 or 3, where 2 means int_mant is wrong by exactly 0.5
689
690    assert(rm != FPRounding_TIEAWAY);
691
692    // Flush to zero:
693    if ((mode & FPLIB_FZ) && exp < 1) {
694        *flags |= FPLIB_UFC;
695        return fp64_zero(sgn);
696    }
697
698    // The bottom 11 bits of mnt are orred together:
699    mnt = (uint64_t)1 << 54 | mnt >> 10 | ((mnt & 0x3ff) != 0);
700
701    if (exp > 0) {
702        biased_exp = exp;
703        int_mant = mnt >> 2;
704        error = mnt & 3;
705    } else {
706        biased_exp = 0;
707        int_mant = lsr64(mnt, 3 - exp);
708        error = (lsr64(mnt, 1 - exp) & 3) | !!(mnt & (lsl64(1, 1 - exp) - 1));
709    }
710
711    if (!biased_exp && error) { // xx should also check fpscr_val<11>
712        *flags |= FPLIB_UFC;
713    }
714
715    // Round up:
716    if ((rm == FPLIB_RN && (error == 3 ||
717                            (error == 2 && (int_mant & 1)))) ||
718        (((rm == FPLIB_RP && !sgn) || (rm == FPLIB_RM && sgn)) && error)) {
719        ++int_mant;
720        if (int_mant == (uint64_t)1 << 52) {
721            // Rounded up from denormalized to normalized
722            biased_exp = 1;
723        }
724        if (int_mant == (uint64_t)1 << 53) {
725            // Rounded up to next exponent
726            ++biased_exp;
727            int_mant >>= 1;
728        }
729    }
730
731    // Handle rounding to odd aka Von Neumann rounding:
732    if (error && rm == FPRounding_ODD)
733        int_mant |= 1;
734
735    // Handle overflow:
736    if (biased_exp >= 2047) {
737        *flags |= FPLIB_OFC | FPLIB_IXC;
738        if (rm == FPLIB_RN || (rm == FPLIB_RP && !sgn) ||
739            (rm == FPLIB_RM && sgn)) {
740            return fp64_infinity(sgn);
741        } else {
742            return fp64_max_normal(sgn);
743        }
744    }
745
746    if (error) {
747        *flags |= FPLIB_IXC;
748    }
749
750    return fp64_pack(sgn, biased_exp, int_mant);
751}
752
753static uint64_t
754fp64_round(int sgn, int exp, uint64_t mnt, int mode, int *flags)
755{
756    return fp64_round_(sgn, exp, mnt, mode & 3, mode, flags);
757}
758
759static int
760fp32_compare_eq(uint32_t a, uint32_t b, int mode, int *flags)
761{
762    int a_sgn, a_exp, b_sgn, b_exp;
763    uint32_t a_mnt, b_mnt;
764
765    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
766    fp32_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
767
768    if ((a_exp == 255 && (uint32_t)(a_mnt << 9)) ||
769        (b_exp == 255 && (uint32_t)(b_mnt << 9))) {
770        if ((a_exp == 255 && (uint32_t)(a_mnt << 9) && !(a >> 22 & 1)) ||
771            (b_exp == 255 && (uint32_t)(b_mnt << 9) && !(b >> 22 & 1)))
772            *flags |= FPLIB_IOC;
773        return 0;
774    }
775    return a == b || (!a_mnt && !b_mnt);
776}
777
778static int
779fp32_compare_ge(uint32_t a, uint32_t b, int mode, int *flags)
780{
781    int a_sgn, a_exp, b_sgn, b_exp;
782    uint32_t a_mnt, b_mnt;
783
784    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
785    fp32_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
786
787    if ((a_exp == 255 && (uint32_t)(a_mnt << 9)) ||
788        (b_exp == 255 && (uint32_t)(b_mnt << 9))) {
789        *flags |= FPLIB_IOC;
790        return 0;
791    }
792    if (!a_mnt && !b_mnt)
793        return 1;
794    if (a_sgn != b_sgn)
795        return b_sgn;
796    if (a_exp != b_exp)
797        return a_sgn ^ (a_exp > b_exp);
798    if (a_mnt != b_mnt)
799        return a_sgn ^ (a_mnt > b_mnt);
800    return 1;
801}
802
803static int
804fp32_compare_gt(uint32_t a, uint32_t b, int mode, int *flags)
805{
806    int a_sgn, a_exp, b_sgn, b_exp;
807    uint32_t a_mnt, b_mnt;
808
809    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
810    fp32_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
811
812    if ((a_exp == 255 && (uint32_t)(a_mnt << 9)) ||
813        (b_exp == 255 && (uint32_t)(b_mnt << 9))) {
814        *flags |= FPLIB_IOC;
815        return 0;
816    }
817    if (!a_mnt && !b_mnt)
818        return 0;
819    if (a_sgn != b_sgn)
820        return b_sgn;
821    if (a_exp != b_exp)
822        return a_sgn ^ (a_exp > b_exp);
823    if (a_mnt != b_mnt)
824        return a_sgn ^ (a_mnt > b_mnt);
825    return 0;
826}
827
828static int
829fp64_compare_eq(uint64_t a, uint64_t b, int mode, int *flags)
830{
831    int a_sgn, a_exp, b_sgn, b_exp;
832    uint64_t a_mnt, b_mnt;
833
834    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
835    fp64_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
836
837    if ((a_exp == 2047 && (uint64_t)(a_mnt << 12)) ||
838        (b_exp == 2047 && (uint64_t)(b_mnt << 12))) {
839        if ((a_exp == 2047 && (uint64_t)(a_mnt << 12) && !(a >> 51 & 1)) ||
840            (b_exp == 2047 && (uint64_t)(b_mnt << 12) && !(b >> 51 & 1)))
841            *flags |= FPLIB_IOC;
842        return 0;
843    }
844    return a == b || (!a_mnt && !b_mnt);
845}
846
847static int
848fp64_compare_ge(uint64_t a, uint64_t b, int mode, int *flags)
849{
850    int a_sgn, a_exp, b_sgn, b_exp;
851    uint64_t a_mnt, b_mnt;
852
853    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
854    fp64_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
855
856    if ((a_exp == 2047 && (uint64_t)(a_mnt << 12)) ||
857        (b_exp == 2047 && (uint64_t)(b_mnt << 12))) {
858        *flags |= FPLIB_IOC;
859        return 0;
860    }
861    if (!a_mnt && !b_mnt)
862        return 1;
863    if (a_sgn != b_sgn)
864        return b_sgn;
865    if (a_exp != b_exp)
866        return a_sgn ^ (a_exp > b_exp);
867    if (a_mnt != b_mnt)
868        return a_sgn ^ (a_mnt > b_mnt);
869    return 1;
870}
871
872static int
873fp64_compare_gt(uint64_t a, uint64_t b, int mode, int *flags)
874{
875    int a_sgn, a_exp, b_sgn, b_exp;
876    uint64_t a_mnt, b_mnt;
877
878    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
879    fp64_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
880
881    if ((a_exp == 2047 && (uint64_t)(a_mnt << 12)) ||
882        (b_exp == 2047 && (uint64_t)(b_mnt << 12))) {
883        *flags |= FPLIB_IOC;
884        return 0;
885    }
886    if (!a_mnt && !b_mnt)
887        return 0;
888    if (a_sgn != b_sgn)
889        return b_sgn;
890    if (a_exp != b_exp)
891        return a_sgn ^ (a_exp > b_exp);
892    if (a_mnt != b_mnt)
893        return a_sgn ^ (a_mnt > b_mnt);
894    return 0;
895}
896
897static uint32_t
898fp32_add(uint32_t a, uint32_t b, int neg, int mode, int *flags)
899{
900    int a_sgn, a_exp, b_sgn, b_exp, x_sgn, x_exp;
901    uint32_t a_mnt, b_mnt, x, x_mnt;
902
903    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
904    fp32_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
905
906    if ((x = fp32_process_NaNs(a, b, mode, flags))) {
907        return x;
908    }
909
910    b_sgn ^= neg;
911
912    // Handle infinities and zeroes:
913    if (a_exp == 255 && b_exp == 255 && a_sgn != b_sgn) {
914        *flags |= FPLIB_IOC;
915        return fp32_defaultNaN();
916    } else if (a_exp == 255) {
917        return fp32_infinity(a_sgn);
918    } else if (b_exp == 255) {
919        return fp32_infinity(b_sgn);
920    } else if (!a_mnt && !b_mnt && a_sgn == b_sgn) {
921        return fp32_zero(a_sgn);
922    }
923
924    a_mnt <<= 3;
925    b_mnt <<= 3;
926    if (a_exp >= b_exp) {
927        b_mnt = (lsr32(b_mnt, a_exp - b_exp) |
928                 !!(b_mnt & (lsl32(1, a_exp - b_exp) - 1)));
929        b_exp = a_exp;
930    } else {
931        a_mnt = (lsr32(a_mnt, b_exp - a_exp) |
932                 !!(a_mnt & (lsl32(1, b_exp - a_exp) - 1)));
933        a_exp = b_exp;
934    }
935    x_sgn = a_sgn;
936    x_exp = a_exp;
937    if (a_sgn == b_sgn) {
938        x_mnt = a_mnt + b_mnt;
939    } else if (a_mnt >= b_mnt) {
940        x_mnt = a_mnt - b_mnt;
941    } else {
942        x_sgn ^= 1;
943        x_mnt = b_mnt - a_mnt;
944    }
945
946    if (!x_mnt) {
947        // Sign of exact zero result depends on rounding mode
948        return fp32_zero((mode & 3) == 2);
949    }
950
951    x_mnt = fp32_normalise(x_mnt, &x_exp);
952
953    return fp32_round(x_sgn, x_exp + 5, x_mnt << 1, mode, flags);
954}
955
956static uint64_t
957fp64_add(uint64_t a, uint64_t b, int neg, int mode, int *flags)
958{
959    int a_sgn, a_exp, b_sgn, b_exp, x_sgn, x_exp;
960    uint64_t a_mnt, b_mnt, x, x_mnt;
961
962    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
963    fp64_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
964
965    if ((x = fp64_process_NaNs(a, b, mode, flags))) {
966        return x;
967    }
968
969    b_sgn ^= neg;
970
971    // Handle infinities and zeroes:
972    if (a_exp == 2047 && b_exp == 2047 && a_sgn != b_sgn) {
973        *flags |= FPLIB_IOC;
974        return fp64_defaultNaN();
975    } else if (a_exp == 2047) {
976        return fp64_infinity(a_sgn);
977    } else if (b_exp == 2047) {
978        return fp64_infinity(b_sgn);
979    } else if (!a_mnt && !b_mnt && a_sgn == b_sgn) {
980        return fp64_zero(a_sgn);
981    }
982
983    a_mnt <<= 3;
984    b_mnt <<= 3;
985    if (a_exp >= b_exp) {
986        b_mnt = (lsr64(b_mnt, a_exp - b_exp) |
987                 !!(b_mnt & (lsl64(1, a_exp - b_exp) - 1)));
988        b_exp = a_exp;
989    } else {
990        a_mnt = (lsr64(a_mnt, b_exp - a_exp) |
991                 !!(a_mnt & (lsl64(1, b_exp - a_exp) - 1)));
992        a_exp = b_exp;
993    }
994    x_sgn = a_sgn;
995    x_exp = a_exp;
996    if (a_sgn == b_sgn) {
997        x_mnt = a_mnt + b_mnt;
998    } else if (a_mnt >= b_mnt) {
999        x_mnt = a_mnt - b_mnt;
1000    } else {
1001        x_sgn ^= 1;
1002        x_mnt = b_mnt - a_mnt;
1003    }
1004
1005    if (!x_mnt) {
1006        // Sign of exact zero result depends on rounding mode
1007        return fp64_zero((mode & 3) == 2);
1008    }
1009
1010    x_mnt = fp64_normalise(x_mnt, &x_exp);
1011
1012    return fp64_round(x_sgn, x_exp + 8, x_mnt << 1, mode, flags);
1013}
1014
1015static uint32_t
1016fp32_mul(uint32_t a, uint32_t b, int mode, int *flags)
1017{
1018    int a_sgn, a_exp, b_sgn, b_exp, x_sgn, x_exp;
1019    uint32_t a_mnt, b_mnt, x;
1020    uint64_t x_mnt;
1021
1022    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1023    fp32_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1024
1025    if ((x = fp32_process_NaNs(a, b, mode, flags))) {
1026        return x;
1027    }
1028
1029    // Handle infinities and zeroes:
1030    if ((a_exp == 255 && !b_mnt) || (b_exp == 255 && !a_mnt)) {
1031        *flags |= FPLIB_IOC;
1032        return fp32_defaultNaN();
1033    } else if (a_exp == 255 || b_exp == 255) {
1034        return fp32_infinity(a_sgn ^ b_sgn);
1035    } else if (!a_mnt || !b_mnt) {
1036        return fp32_zero(a_sgn ^ b_sgn);
1037    }
1038
1039    // Multiply and normalise:
1040    x_sgn = a_sgn ^ b_sgn;
1041    x_exp = a_exp + b_exp - 110;
1042    x_mnt = (uint64_t)a_mnt * b_mnt;
1043    x_mnt = fp64_normalise(x_mnt, &x_exp);
1044
1045    // Convert to 32 bits, collapsing error into bottom bit:
1046    x_mnt = lsr64(x_mnt, 31) | !!lsl64(x_mnt, 33);
1047
1048    return fp32_round(x_sgn, x_exp, x_mnt, mode, flags);
1049}
1050
1051static uint64_t
1052fp64_mul(uint64_t a, uint64_t b, int mode, int *flags)
1053{
1054    int a_sgn, a_exp, b_sgn, b_exp, x_sgn, x_exp;
1055    uint64_t a_mnt, b_mnt, x;
1056    uint64_t x0_mnt, x1_mnt;
1057
1058    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1059    fp64_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1060
1061    if ((x = fp64_process_NaNs(a, b, mode, flags))) {
1062        return x;
1063    }
1064
1065    // Handle infinities and zeroes:
1066    if ((a_exp == 2047 && !b_mnt) || (b_exp == 2047 && !a_mnt)) {
1067        *flags |= FPLIB_IOC;
1068        return fp64_defaultNaN();
1069    } else if (a_exp == 2047 || b_exp == 2047) {
1070        return fp64_infinity(a_sgn ^ b_sgn);
1071    } else if (!a_mnt || !b_mnt) {
1072        return fp64_zero(a_sgn ^ b_sgn);
1073    }
1074
1075    // Multiply and normalise:
1076    x_sgn = a_sgn ^ b_sgn;
1077    x_exp = a_exp + b_exp - 1000;
1078    mul62x62(&x0_mnt, &x1_mnt, a_mnt, b_mnt);
1079    fp128_normalise(&x0_mnt, &x1_mnt, &x_exp);
1080
1081    // Convert to 64 bits, collapsing error into bottom bit:
1082    x0_mnt = x1_mnt << 1 | !!x0_mnt;
1083
1084    return fp64_round(x_sgn, x_exp, x0_mnt, mode, flags);
1085}
1086
1087static uint32_t
1088fp32_muladd(uint32_t a, uint32_t b, uint32_t c, int scale,
1089            int mode, int *flags)
1090{
1091    int a_sgn, a_exp, b_sgn, b_exp, c_sgn, c_exp, x_sgn, x_exp, y_sgn, y_exp;
1092    uint32_t a_mnt, b_mnt, c_mnt, x;
1093    uint64_t x_mnt, y_mnt;
1094
1095    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1096    fp32_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1097    fp32_unpack(&c_sgn, &c_exp, &c_mnt, c, mode, flags);
1098
1099    x = fp32_process_NaNs3(a, b, c, mode, flags);
1100
1101    // Quiet NaN added to product of zero and infinity:
1102    if (a_exp == 255 && (a_mnt >> 22 & 1) &&
1103        ((!b_mnt && c_exp == 255 && !(uint32_t)(c_mnt << 9)) ||
1104         (!c_mnt && b_exp == 255 && !(uint32_t)(b_mnt << 9)))) {
1105        x = fp32_defaultNaN();
1106        *flags |= FPLIB_IOC;
1107    }
1108
1109    if (x) {
1110        return x;
1111    }
1112
1113    // Handle infinities and zeroes:
1114    if ((b_exp == 255 && !c_mnt) ||
1115        (c_exp == 255 && !b_mnt) ||
1116        (a_exp == 255 && (b_exp == 255 || c_exp == 255) &&
1117         (a_sgn != (b_sgn ^ c_sgn)))) {
1118        *flags |= FPLIB_IOC;
1119        return fp32_defaultNaN();
1120    }
1121    if (a_exp == 255)
1122        return fp32_infinity(a_sgn);
1123    if (b_exp == 255 || c_exp == 255)
1124        return fp32_infinity(b_sgn ^ c_sgn);
1125    if (!a_mnt && (!b_mnt || !c_mnt) && a_sgn == (b_sgn ^ c_sgn))
1126        return fp32_zero(a_sgn);
1127
1128    x_sgn = a_sgn;
1129    x_exp = a_exp + 13;
1130    x_mnt = (uint64_t)a_mnt << 27;
1131
1132    // Multiply:
1133    y_sgn = b_sgn ^ c_sgn;
1134    y_exp = b_exp + c_exp - 113;
1135    y_mnt = (uint64_t)b_mnt * c_mnt << 3;
1136    if (!y_mnt) {
1137        y_exp = x_exp;
1138    }
1139
1140    // Add:
1141    if (x_exp >= y_exp) {
1142        y_mnt = (lsr64(y_mnt, x_exp - y_exp) |
1143                 !!(y_mnt & (lsl64(1, x_exp - y_exp) - 1)));
1144        y_exp = x_exp;
1145    } else {
1146        x_mnt = (lsr64(x_mnt, y_exp - x_exp) |
1147                 !!(x_mnt & (lsl64(1, y_exp - x_exp) - 1)));
1148        x_exp = y_exp;
1149    }
1150    if (x_sgn == y_sgn) {
1151        x_mnt = x_mnt + y_mnt;
1152    } else if (x_mnt >= y_mnt) {
1153        x_mnt = x_mnt - y_mnt;
1154    } else {
1155        x_sgn ^= 1;
1156        x_mnt = y_mnt - x_mnt;
1157    }
1158
1159    if (!x_mnt) {
1160        // Sign of exact zero result depends on rounding mode
1161        return fp32_zero((mode & 3) == 2);
1162    }
1163
1164    // Normalise and convert to 32 bits, collapsing error into bottom bit:
1165    x_mnt = fp64_normalise(x_mnt, &x_exp);
1166    x_mnt = x_mnt >> 31 | !!(uint32_t)(x_mnt << 1);
1167
1168    return fp32_round(x_sgn, x_exp + scale, x_mnt, mode, flags);
1169}
1170
1171static uint64_t
1172fp64_muladd(uint64_t a, uint64_t b, uint64_t c, int scale,
1173            int mode, int *flags)
1174{
1175    int a_sgn, a_exp, b_sgn, b_exp, c_sgn, c_exp, x_sgn, x_exp, y_sgn, y_exp;
1176    uint64_t a_mnt, b_mnt, c_mnt, x;
1177    uint64_t x0_mnt, x1_mnt, y0_mnt, y1_mnt;
1178
1179    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1180    fp64_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1181    fp64_unpack(&c_sgn, &c_exp, &c_mnt, c, mode, flags);
1182
1183    x = fp64_process_NaNs3(a, b, c, mode, flags);
1184
1185    // Quiet NaN added to product of zero and infinity:
1186    if (a_exp == 2047 && (a_mnt >> 51 & 1) &&
1187        ((!b_mnt && c_exp == 2047 && !(uint64_t)(c_mnt << 12)) ||
1188         (!c_mnt && b_exp == 2047 && !(uint64_t)(b_mnt << 12)))) {
1189        x = fp64_defaultNaN();
1190        *flags |= FPLIB_IOC;
1191    }
1192
1193    if (x) {
1194        return x;
1195    }
1196
1197    // Handle infinities and zeroes:
1198    if ((b_exp == 2047 && !c_mnt) ||
1199        (c_exp == 2047 && !b_mnt) ||
1200        (a_exp == 2047 && (b_exp == 2047 || c_exp == 2047) &&
1201         (a_sgn != (b_sgn ^ c_sgn)))) {
1202        *flags |= FPLIB_IOC;
1203        return fp64_defaultNaN();
1204    }
1205    if (a_exp == 2047)
1206        return fp64_infinity(a_sgn);
1207    if (b_exp == 2047 || c_exp == 2047)
1208        return fp64_infinity(b_sgn ^ c_sgn);
1209    if (!a_mnt && (!b_mnt || !c_mnt) && a_sgn == (b_sgn ^ c_sgn))
1210        return fp64_zero(a_sgn);
1211
1212    x_sgn = a_sgn;
1213    x_exp = a_exp + 11;
1214    x0_mnt = 0;
1215    x1_mnt = a_mnt;
1216
1217    // Multiply:
1218    y_sgn = b_sgn ^ c_sgn;
1219    y_exp = b_exp + c_exp - 1003;
1220    mul62x62(&y0_mnt, &y1_mnt, b_mnt, c_mnt << 3);
1221    if (!y0_mnt && !y1_mnt) {
1222        y_exp = x_exp;
1223    }
1224
1225    // Add:
1226    if (x_exp >= y_exp) {
1227        uint64_t t0, t1;
1228        lsl128(&t0, &t1, y0_mnt, y1_mnt,
1229               x_exp - y_exp < 128 ? 128 - (x_exp - y_exp) : 0);
1230        lsr128(&y0_mnt, &y1_mnt, y0_mnt, y1_mnt, x_exp - y_exp);
1231        y0_mnt |= !!(t0 | t1);
1232        y_exp = x_exp;
1233    } else {
1234        uint64_t t0, t1;
1235        lsl128(&t0, &t1, x0_mnt, x1_mnt,
1236               y_exp - x_exp < 128 ? 128 - (y_exp - x_exp) : 0);
1237        lsr128(&x0_mnt, &x1_mnt, x0_mnt, x1_mnt, y_exp - x_exp);
1238        x0_mnt |= !!(t0 | t1);
1239        x_exp = y_exp;
1240    }
1241    if (x_sgn == y_sgn) {
1242        add128(&x0_mnt, &x1_mnt, x0_mnt, x1_mnt, y0_mnt, y1_mnt);
1243    } else if (cmp128(x0_mnt, x1_mnt, y0_mnt, y1_mnt) >= 0) {
1244        sub128(&x0_mnt, &x1_mnt, x0_mnt, x1_mnt, y0_mnt, y1_mnt);
1245    } else {
1246        x_sgn ^= 1;
1247        sub128(&x0_mnt, &x1_mnt, y0_mnt, y1_mnt, x0_mnt, x1_mnt);
1248    }
1249
1250    if (!x0_mnt && !x1_mnt) {
1251        // Sign of exact zero result depends on rounding mode
1252        return fp64_zero((mode & 3) == 2);
1253    }
1254
1255    // Normalise and convert to 64 bits, collapsing error into bottom bit:
1256    fp128_normalise(&x0_mnt, &x1_mnt, &x_exp);
1257    x0_mnt = x1_mnt << 1 | !!x0_mnt;
1258
1259    return fp64_round(x_sgn, x_exp + scale, x0_mnt, mode, flags);
1260}
1261
1262static uint32_t
1263fp32_div(uint32_t a, uint32_t b, int mode, int *flags)
1264{
1265    int a_sgn, a_exp, b_sgn, b_exp, x_sgn, x_exp;
1266    uint32_t a_mnt, b_mnt, x;
1267    uint64_t x_mnt;
1268
1269    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1270    fp32_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1271
1272    if ((x = fp32_process_NaNs(a, b, mode, flags)))
1273        return x;
1274
1275    // Handle infinities and zeroes:
1276    if ((a_exp == 255 && b_exp == 255) || (!a_mnt && !b_mnt)) {
1277        *flags |= FPLIB_IOC;
1278        return fp32_defaultNaN();
1279    }
1280    if (a_exp == 255 || !b_mnt) {
1281        if (a_exp != 255)
1282            *flags |= FPLIB_DZC;
1283        return fp32_infinity(a_sgn ^ b_sgn);
1284    }
1285    if (!a_mnt || b_exp == 255)
1286        return fp32_zero(a_sgn ^ b_sgn);
1287
1288    // Divide, setting bottom bit if inexact:
1289    a_mnt = fp32_normalise(a_mnt, &a_exp);
1290    x_sgn = a_sgn ^ b_sgn;
1291    x_exp = a_exp - b_exp + 172;
1292    x_mnt = ((uint64_t)a_mnt << 18) / b_mnt;
1293    x_mnt |= (x_mnt * b_mnt != (uint64_t)a_mnt << 18);
1294
1295    // Normalise and convert to 32 bits, collapsing error into bottom bit:
1296    x_mnt = fp64_normalise(x_mnt, &x_exp);
1297    x_mnt = x_mnt >> 31 | !!(uint32_t)(x_mnt << 1);
1298
1299    return fp32_round(x_sgn, x_exp, x_mnt, mode, flags);
1300}
1301
1302static uint64_t
1303fp64_div(uint64_t a, uint64_t b, int mode, int *flags)
1304{
1305    int a_sgn, a_exp, b_sgn, b_exp, x_sgn, x_exp, c;
1306    uint64_t a_mnt, b_mnt, x, x_mnt, x0_mnt, x1_mnt;
1307
1308    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1309    fp64_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1310
1311    if ((x = fp64_process_NaNs(a, b, mode, flags)))
1312        return x;
1313
1314    // Handle infinities and zeroes:
1315    if ((a_exp == 2047 && b_exp == 2047) || (!a_mnt && !b_mnt)) {
1316        *flags |= FPLIB_IOC;
1317        return fp64_defaultNaN();
1318    }
1319    if (a_exp == 2047 || !b_mnt) {
1320        if (a_exp != 2047)
1321            *flags |= FPLIB_DZC;
1322        return fp64_infinity(a_sgn ^ b_sgn);
1323    }
1324    if (!a_mnt || b_exp == 2047)
1325        return fp64_zero(a_sgn ^ b_sgn);
1326
1327    // Find reciprocal of divisor with Newton-Raphson:
1328    a_mnt = fp64_normalise(a_mnt, &a_exp);
1329    b_mnt = fp64_normalise(b_mnt, &b_exp);
1330    x_mnt = ~(uint64_t)0 / (b_mnt >> 31);
1331    mul64x32(&x0_mnt, &x1_mnt, b_mnt, x_mnt);
1332    sub128(&x0_mnt, &x1_mnt, 0, (uint64_t)1 << 32, x0_mnt, x1_mnt);
1333    lsr128(&x0_mnt, &x1_mnt, x0_mnt, x1_mnt, 32);
1334    mul64x32(&x0_mnt, &x1_mnt, x0_mnt, x_mnt);
1335    lsr128(&x0_mnt, &x1_mnt, x0_mnt, x1_mnt, 33);
1336
1337    // Multiply by dividend:
1338    x_sgn = a_sgn ^ b_sgn;
1339    x_exp = a_exp - b_exp + 1031;
1340    mul62x62(&x0_mnt, &x1_mnt, x0_mnt, a_mnt >> 2); // xx 62x62 is enough
1341    lsr128(&x0_mnt, &x1_mnt, x0_mnt, x1_mnt, 4);
1342    x_mnt = x1_mnt;
1343
1344    // This is an underestimate, so try adding one:
1345    mul62x62(&x0_mnt, &x1_mnt, b_mnt >> 2, x_mnt + 1); // xx 62x62 is enough
1346    c = cmp128(x0_mnt, x1_mnt, 0, a_mnt >> 11);
1347    if (c <= 0) {
1348        ++x_mnt;
1349    }
1350
1351    x_mnt = fp64_normalise(x_mnt, &x_exp);
1352
1353    return fp64_round(x_sgn, x_exp, x_mnt << 1 | !!c, mode, flags);
1354}
1355
1356static void
1357set_fpscr0(FPSCR &fpscr, int flags)
1358{
1359    if (flags & FPLIB_IDC) {
1360        fpscr.idc = 1;
1361    }
1362    if (flags & FPLIB_IOC) {
1363        fpscr.ioc = 1;
1364    }
1365    if (flags & FPLIB_DZC) {
1366        fpscr.dzc = 1;
1367    }
1368    if (flags & FPLIB_OFC) {
1369        fpscr.ofc = 1;
1370    }
1371    if (flags & FPLIB_UFC) {
1372        fpscr.ufc = 1;
1373    }
1374    if (flags & FPLIB_IXC) {
1375        fpscr.ixc = 1;
1376    }
1377}
1378
1379static uint32_t
1380fp32_sqrt(uint32_t a, int mode, int *flags)
1381{
1382    int a_sgn, a_exp, x_sgn, x_exp;
1383    uint32_t a_mnt, x, x_mnt;
1384    uint64_t t0, t1;
1385
1386    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1387
1388    // Handle NaNs:
1389    if (a_exp == 255 && (uint32_t)(a_mnt << 9))
1390        return fp32_process_NaN(a, mode, flags);
1391
1392    // Handle infinities and zeroes:
1393    if (!a_mnt) {
1394        return fp32_zero(a_sgn);
1395    }
1396    if (a_exp == 255 && !a_sgn) {
1397        return fp32_infinity(a_sgn);
1398    }
1399    if (a_sgn) {
1400        *flags |= FPLIB_IOC;
1401        return fp32_defaultNaN();
1402    }
1403
1404    a_mnt = fp32_normalise(a_mnt, &a_exp);
1405    if (!(a_exp & 1)) {
1406        ++a_exp;
1407        a_mnt >>= 1;
1408    }
1409
1410    // x = (a * 3 + 5) / 8
1411    x = (a_mnt >> 2) + (a_mnt >> 3) + (5 << 28);
1412
1413    // x = (a / x + x) / 2; // 16-bit accuracy
1414    x = (a_mnt / (x >> 15) + (x >> 16)) << 15;
1415
1416    // x = (a / x + x) / 2; // 16-bit accuracy
1417    x = (a_mnt / (x >> 15) + (x >> 16)) << 15;
1418
1419    // x = (a / x + x) / 2; // 32-bit accuracy
1420    x = ((((uint64_t)a_mnt << 32) / x) >> 2) + (x >> 1);
1421
1422    x_sgn = 0;
1423    x_exp = (a_exp + 147) >> 1;
1424    x_mnt = ((x - (1 << 5)) >> 6) + 1;
1425    t1 = (uint64_t)x_mnt * x_mnt;
1426    t0 = (uint64_t)a_mnt << 19;
1427    if (t1 > t0) {
1428        --x_mnt;
1429    }
1430
1431    x_mnt = fp32_normalise(x_mnt, &x_exp);
1432
1433    return fp32_round(x_sgn, x_exp, x_mnt << 1 | (t1 != t0), mode, flags);
1434}
1435
1436static uint64_t
1437fp64_sqrt(uint64_t a, int mode, int *flags)
1438{
1439    int a_sgn, a_exp, x_sgn, x_exp, c;
1440    uint64_t a_mnt, x_mnt, r, x0, x1;
1441    uint32_t x;
1442
1443    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1444
1445    // Handle NaNs:
1446    if (a_exp == 2047 && (uint64_t)(a_mnt << 12)) {
1447        return fp64_process_NaN(a, mode, flags);
1448    }
1449
1450    // Handle infinities and zeroes:
1451    if (!a_mnt)
1452        return fp64_zero(a_sgn);
1453    if (a_exp == 2047 && !a_sgn)
1454        return fp64_infinity(a_sgn);
1455    if (a_sgn) {
1456        *flags |= FPLIB_IOC;
1457        return fp64_defaultNaN();
1458    }
1459
1460    a_mnt = fp64_normalise(a_mnt, &a_exp);
1461    if (a_exp & 1) {
1462        ++a_exp;
1463        a_mnt >>= 1;
1464    }
1465
1466    // x = (a * 3 + 5) / 8
1467    x = (a_mnt >> 34) + (a_mnt >> 35) + (5 << 28);
1468
1469    // x = (a / x + x) / 2; // 16-bit accuracy
1470    x = ((a_mnt >> 32) / (x >> 15) + (x >> 16)) << 15;
1471
1472    // x = (a / x + x) / 2; // 16-bit accuracy
1473    x = ((a_mnt >> 32) / (x >> 15) + (x >> 16)) << 15;
1474
1475    // x = (a / x + x) / 2; // 32-bit accuracy
1476    x = ((a_mnt / x) >> 2) + (x >> 1);
1477
1478    // r = 1 / x; // 32-bit accuracy
1479    r = ((uint64_t)1 << 62) / x;
1480
1481    // r = r * (2 - x * r); // 64-bit accuracy
1482    mul64x32(&x0, &x1, -(uint64_t)x * r << 1, r);
1483    lsr128(&x0, &x1, x0, x1, 31);
1484
1485    // x = (x + a * r) / 2; // 64-bit accuracy
1486    mul62x62(&x0, &x1, a_mnt >> 10, x0 >> 2);
1487    lsl128(&x0, &x1, x0, x1, 5);
1488    lsr128(&x0, &x1, x0, x1, 56);
1489
1490    x0 = ((uint64_t)x << 31) + (x0 >> 1);
1491
1492    x_sgn = 0;
1493    x_exp = (a_exp + 1053) >> 1;
1494    x_mnt = x0;
1495    x_mnt = ((x_mnt - (1 << 8)) >> 9) + 1;
1496    mul62x62(&x0, &x1, x_mnt, x_mnt);
1497    lsl128(&x0, &x1, x0, x1, 19);
1498    c = cmp128(x0, x1, 0, a_mnt);
1499    if (c > 0)
1500        --x_mnt;
1501
1502    x_mnt = fp64_normalise(x_mnt, &x_exp);
1503
1504    return fp64_round(x_sgn, x_exp, x_mnt << 1 | !!c, mode, flags);
1505}
1506
1507static int
1508modeConv(FPSCR fpscr)
1509{
1510    return (((int) fpscr) >> 22) & 0xF;
1511}
1512
1513static void
1514set_fpscr(FPSCR &fpscr, int flags)
1515{
1516    // translate back to FPSCR
1517    bool underflow = false;
1518    if (flags & FPLIB_IDC) {
1519        fpscr.idc = 1;
1520    }
1521    if (flags & FPLIB_IOC) {
1522        fpscr.ioc = 1;
1523    }
1524    if (flags & FPLIB_DZC) {
1525        fpscr.dzc = 1;
1526    }
1527    if (flags & FPLIB_OFC) {
1528        fpscr.ofc = 1;
1529    }
1530    if (flags & FPLIB_UFC) {
1531        underflow = true; //xx Why is this required?
1532        fpscr.ufc = 1;
1533    }
1534    if ((flags & FPLIB_IXC) && !(underflow && fpscr.fz)) {
1535        fpscr.ixc = 1;
1536    }
1537}
1538
1539template <>
1540bool
1541fplibCompareEQ(uint32_t a, uint32_t b, FPSCR &fpscr)
1542{
1543    int flags = 0;
1544    int x = fp32_compare_eq(a, b, modeConv(fpscr), &flags);
1545    set_fpscr(fpscr, flags);
1546    return x;
1547}
1548
1549template <>
1550bool
1551fplibCompareGE(uint32_t a, uint32_t b, FPSCR &fpscr)
1552{
1553    int flags = 0;
1554    int x = fp32_compare_ge(a, b, modeConv(fpscr), &flags);
1555    set_fpscr(fpscr, flags);
1556    return x;
1557}
1558
1559template <>
1560bool
1561fplibCompareGT(uint32_t a, uint32_t b, FPSCR &fpscr)
1562{
1563    int flags = 0;
1564    int x = fp32_compare_gt(a, b, modeConv(fpscr), &flags);
1565    set_fpscr(fpscr, flags);
1566    return x;
1567}
1568
1569template <>
1570bool
1571fplibCompareEQ(uint64_t a, uint64_t b, FPSCR &fpscr)
1572{
1573    int flags = 0;
1574    int x = fp64_compare_eq(a, b, modeConv(fpscr), &flags);
1575    set_fpscr(fpscr, flags);
1576    return x;
1577}
1578
1579template <>
1580bool
1581fplibCompareGE(uint64_t a, uint64_t b, FPSCR &fpscr)
1582{
1583    int flags = 0;
1584    int x = fp64_compare_ge(a, b, modeConv(fpscr), &flags);
1585    set_fpscr(fpscr, flags);
1586    return x;
1587}
1588
1589template <>
1590bool
1591fplibCompareGT(uint64_t a, uint64_t b, FPSCR &fpscr)
1592{
1593    int flags = 0;
1594    int x = fp64_compare_gt(a, b, modeConv(fpscr), &flags);
1595    set_fpscr(fpscr, flags);
1596    return x;
1597}
1598
1599template <>
1600uint32_t
1601fplibAbs(uint32_t op)
1602{
1603    return op & ~((uint32_t)1 << 31);
1604}
1605
1606template <>
1607uint64_t
1608fplibAbs(uint64_t op)
1609{
1610    return op & ~((uint64_t)1 << 63);
1611}
1612
1613template <>
1614uint32_t
1615fplibAdd(uint32_t op1, uint32_t op2, FPSCR &fpscr)
1616{
1617    int flags = 0;
1618    uint32_t result = fp32_add(op1, op2, 0, modeConv(fpscr), &flags);
1619    set_fpscr0(fpscr, flags);
1620    return result;
1621}
1622
1623template <>
1624uint64_t
1625fplibAdd(uint64_t op1, uint64_t op2, FPSCR &fpscr)
1626{
1627    int flags = 0;
1628    uint64_t result = fp64_add(op1, op2, 0, modeConv(fpscr), &flags);
1629    set_fpscr0(fpscr, flags);
1630    return result;
1631}
1632
1633template <>
1634int
1635fplibCompare(uint32_t op1, uint32_t op2, bool signal_nans, FPSCR &fpscr)
1636{
1637    int mode = modeConv(fpscr);
1638    int flags = 0;
1639    int sgn1, exp1, sgn2, exp2, result;
1640    uint32_t mnt1, mnt2;
1641
1642    fp32_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
1643    fp32_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
1644
1645    if ((exp1 == 255 && (uint32_t)(mnt1 << 9)) ||
1646        (exp2 == 255 && (uint32_t)(mnt2 << 9))) {
1647        result = 3;
1648        if ((exp1 == 255 && (uint32_t)(mnt1 << 9) && !(mnt1 >> 22 & 1)) ||
1649            (exp2 == 255 && (uint32_t)(mnt2 << 9) && !(mnt2 >> 22 & 1)) ||
1650            signal_nans)
1651            flags |= FPLIB_IOC;
1652    } else {
1653        if (op1 == op2 || (!mnt1 && !mnt2)) {
1654            result = 6;
1655        } else if (sgn1 != sgn2) {
1656            result = sgn1 ? 8 : 2;
1657        } else if (exp1 != exp2) {
1658            result = sgn1 ^ (exp1 < exp2) ? 8 : 2;
1659        } else {
1660            result = sgn1 ^ (mnt1 < mnt2) ? 8 : 2;
1661        }
1662    }
1663
1664    set_fpscr0(fpscr, flags);
1665
1666    return result;
1667}
1668
1669template <>
1670int
1671fplibCompare(uint64_t op1, uint64_t op2, bool signal_nans, FPSCR &fpscr)
1672{
1673    int mode = modeConv(fpscr);
1674    int flags = 0;
1675    int sgn1, exp1, sgn2, exp2, result;
1676    uint64_t mnt1, mnt2;
1677
1678    fp64_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
1679    fp64_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
1680
1681    if ((exp1 == 2047 && (uint64_t)(mnt1 << 12)) ||
1682        (exp2 == 2047 && (uint64_t)(mnt2 << 12))) {
1683        result = 3;
1684        if ((exp1 == 2047 && (uint64_t)(mnt1 << 12) && !(mnt1 >> 51 & 1)) ||
1685            (exp2 == 2047 && (uint64_t)(mnt2 << 12) && !(mnt2 >> 51 & 1)) ||
1686            signal_nans)
1687            flags |= FPLIB_IOC;
1688    } else {
1689        if (op1 == op2 || (!mnt1 && !mnt2)) {
1690            result = 6;
1691        } else if (sgn1 != sgn2) {
1692            result = sgn1 ? 8 : 2;
1693        } else if (exp1 != exp2) {
1694            result = sgn1 ^ (exp1 < exp2) ? 8 : 2;
1695        } else {
1696            result = sgn1 ^ (mnt1 < mnt2) ? 8 : 2;
1697        }
1698    }
1699
1700    set_fpscr0(fpscr, flags);
1701
1702    return result;
1703}
1704
1705static uint16_t
1706fp16_FPConvertNaN_32(uint32_t op)
1707{
1708    return fp16_pack(op >> 31, 31, (uint16_t)1 << 9 | op >> 13);
1709}
1710
1711static uint16_t
1712fp16_FPConvertNaN_64(uint64_t op)
1713{
1714    return fp16_pack(op >> 63, 31, (uint16_t)1 << 9 | op >> 42);
1715}
1716
1717static uint32_t
1718fp32_FPConvertNaN_16(uint16_t op)
1719{
1720    return fp32_pack(op >> 15, 255, (uint32_t)1 << 22 | (uint32_t)op << 13);
1721}
1722
1723static uint32_t
1724fp32_FPConvertNaN_64(uint64_t op)
1725{
1726    return fp32_pack(op >> 63, 255, (uint32_t)1 << 22 | op >> 29);
1727}
1728
1729static uint64_t
1730fp64_FPConvertNaN_16(uint16_t op)
1731{
1732    return fp64_pack(op >> 15, 2047, (uint64_t)1 << 51 | (uint64_t)op << 42);
1733}
1734
1735static uint64_t
1736fp64_FPConvertNaN_32(uint32_t op)
1737{
1738    return fp64_pack(op >> 31, 2047, (uint64_t)1 << 51 | (uint64_t)op << 29);
1739}
1740
1741static uint32_t
1742fp32_FPOnePointFive(int sgn)
1743{
1744    return fp32_pack(sgn, 127, (uint64_t)1 << 22);
1745}
1746
1747static uint64_t
1748fp64_FPOnePointFive(int sgn)
1749{
1750    return fp64_pack(sgn, 1023, (uint64_t)1 << 51);
1751}
1752
1753static uint32_t
1754fp32_FPThree(int sgn)
1755{
1756    return fp32_pack(sgn, 128, (uint64_t)1 << 22);
1757}
1758
1759static uint64_t
1760fp64_FPThree(int sgn)
1761{
1762    return fp64_pack(sgn, 1024, (uint64_t)1 << 51);
1763}
1764
1765static uint32_t
1766fp32_FPTwo(int sgn)
1767{
1768    return fp32_pack(sgn, 128, 0);
1769}
1770
1771static uint64_t
1772fp64_FPTwo(int sgn)
1773{
1774    return fp64_pack(sgn, 1024, 0);
1775}
1776
1777template <>
1778uint16_t
1779fplibConvert(uint32_t op, FPRounding rounding, FPSCR &fpscr)
1780{
1781    int mode = modeConv(fpscr);
1782    int flags = 0;
1783    int sgn, exp;
1784    uint32_t mnt;
1785    uint16_t result;
1786
1787    // Unpack floating-point operand optionally with flush-to-zero:
1788    fp32_unpack(&sgn, &exp, &mnt, op, mode, &flags);
1789
1790    bool alt_hp = fpscr.ahp;
1791
1792    if (exp == 255 && (uint32_t)(mnt << 9)) {
1793        if (alt_hp) {
1794            result = fp16_zero(sgn);
1795        } else if (fpscr.dn) {
1796            result = fp16_defaultNaN();
1797        } else {
1798            result = fp16_FPConvertNaN_32(op);
1799        }
1800        if (!(mnt >> 22 & 1) || alt_hp) {
1801            flags |= FPLIB_IOC;
1802        }
1803    } else if (exp == 255) {
1804        if (alt_hp) {
1805            result = sgn << 15 | (uint16_t)0x7fff;
1806            flags |= FPLIB_IOC;
1807        } else {
1808            result = fp16_infinity(sgn);
1809        }
1810    } else if (!mnt) {
1811        result = fp16_zero(sgn);
1812    } else {
1813        result = fp16_round_(sgn, exp - 127 + 15,
1814                             mnt >> 7 | !!(uint32_t)(mnt << 25),
1815                             rounding, mode | alt_hp << 4, &flags);
1816    }
1817
1818    set_fpscr0(fpscr, flags);
1819
1820    return result;
1821}
1822
1823template <>
1824uint16_t
1825fplibConvert(uint64_t op, FPRounding rounding, FPSCR &fpscr)
1826{
1827    int mode = modeConv(fpscr);
1828    int flags = 0;
1829    int sgn, exp;
1830    uint64_t mnt;
1831    uint16_t result;
1832
1833    // Unpack floating-point operand optionally with flush-to-zero:
1834    fp64_unpack(&sgn, &exp, &mnt, op, mode, &flags);
1835
1836    bool alt_hp = fpscr.ahp;
1837
1838    if (exp == 2047 && (uint64_t)(mnt << 12)) {
1839        if (alt_hp) {
1840            result = fp16_zero(sgn);
1841        } else if (fpscr.dn) {
1842            result = fp16_defaultNaN();
1843        } else {
1844            result = fp16_FPConvertNaN_64(op);
1845        }
1846        if (!(mnt >> 51 & 1) || alt_hp) {
1847            flags |= FPLIB_IOC;
1848        }
1849    } else if (exp == 2047) {
1850        if (alt_hp) {
1851            result = sgn << 15 | (uint16_t)0x7fff;
1852            flags |= FPLIB_IOC;
1853        } else {
1854            result = fp16_infinity(sgn);
1855        }
1856    } else if (!mnt) {
1857        result = fp16_zero(sgn);
1858    } else {
1859        result = fp16_round_(sgn, exp - 1023 + 15,
1860                             mnt >> 36 | !!(uint64_t)(mnt << 28),
1861                             rounding, mode | alt_hp << 4, &flags);
1862    }
1863
1864    set_fpscr0(fpscr, flags);
1865
1866    return result;
1867}
1868
1869template <>
1870uint32_t
1871fplibConvert(uint16_t op, FPRounding rounding, FPSCR &fpscr)
1872{
1873    int mode = modeConv(fpscr);
1874    int flags = 0;
1875    int sgn, exp;
1876    uint16_t mnt;
1877    uint32_t result;
1878
1879    // Unpack floating-point operand optionally with flush-to-zero:
1880    fp16_unpack(&sgn, &exp, &mnt, op, mode, &flags);
1881
1882    if (exp == 31 && !fpscr.ahp && (uint16_t)(mnt << 6)) {
1883        if (fpscr.dn) {
1884            result = fp32_defaultNaN();
1885        } else {
1886            result = fp32_FPConvertNaN_16(op);
1887        }
1888        if (!(mnt >> 9 & 1)) {
1889            flags |= FPLIB_IOC;
1890        }
1891    } else if (exp == 31 && !fpscr.ahp) {
1892        result = fp32_infinity(sgn);
1893    } else if (!mnt) {
1894        result = fp32_zero(sgn);
1895    } else {
1896        mnt = fp16_normalise(mnt, &exp);
1897        result = fp32_pack(sgn, exp - 15 + 127 + 5, (uint32_t)mnt << 8);
1898    }
1899
1900    set_fpscr0(fpscr, flags);
1901
1902    return result;
1903}
1904
1905template <>
1906uint32_t
1907fplibConvert(uint64_t op, FPRounding rounding, FPSCR &fpscr)
1908{
1909    int mode = modeConv(fpscr);
1910    int flags = 0;
1911    int sgn, exp;
1912    uint64_t mnt;
1913    uint32_t result;
1914
1915    // Unpack floating-point operand optionally with flush-to-zero:
1916    fp64_unpack(&sgn, &exp, &mnt, op, mode, &flags);
1917
1918    if (exp == 2047 && (uint64_t)(mnt << 12)) {
1919        if (fpscr.dn) {
1920            result = fp32_defaultNaN();
1921        } else {
1922            result = fp32_FPConvertNaN_64(op);
1923        }
1924        if (!(mnt >> 51 & 1)) {
1925            flags |= FPLIB_IOC;
1926        }
1927    } else if (exp == 2047) {
1928        result = fp32_infinity(sgn);
1929    } else if (!mnt) {
1930        result = fp32_zero(sgn);
1931    } else {
1932        result = fp32_round_(sgn, exp - 1023 + 127,
1933                             mnt >> 20 | !!(uint64_t)(mnt << 44),
1934                             rounding, mode, &flags);
1935    }
1936
1937    set_fpscr0(fpscr, flags);
1938
1939    return result;
1940}
1941
1942template <>
1943uint64_t
1944fplibConvert(uint16_t op, FPRounding rounding, FPSCR &fpscr)
1945{
1946    int mode = modeConv(fpscr);
1947    int flags = 0;
1948    int sgn, exp;
1949    uint16_t mnt;
1950    uint64_t result;
1951
1952    // Unpack floating-point operand optionally with flush-to-zero:
1953    fp16_unpack(&sgn, &exp, &mnt, op, mode, &flags);
1954
1955    if (exp == 31 && !fpscr.ahp && (uint16_t)(mnt << 6)) {
1956        if (fpscr.dn) {
1957            result = fp64_defaultNaN();
1958        } else {
1959            result = fp64_FPConvertNaN_16(op);
1960        }
1961        if (!(mnt >> 9 & 1)) {
1962            flags |= FPLIB_IOC;
1963        }
1964    } else if (exp == 31 && !fpscr.ahp) {
1965        result = fp64_infinity(sgn);
1966    } else if (!mnt) {
1967        result = fp64_zero(sgn);
1968    } else {
1969        mnt = fp16_normalise(mnt, &exp);
1970        result = fp64_pack(sgn, exp - 15 + 1023 + 5, (uint64_t)mnt << 37);
1971    }
1972
1973    set_fpscr0(fpscr, flags);
1974
1975    return result;
1976}
1977
1978template <>
1979uint64_t
1980fplibConvert(uint32_t op, FPRounding rounding, FPSCR &fpscr)
1981{
1982    int mode = modeConv(fpscr);
1983    int flags = 0;
1984    int sgn, exp;
1985    uint32_t mnt;
1986    uint64_t result;
1987
1988    // Unpack floating-point operand optionally with flush-to-zero:
1989    fp32_unpack(&sgn, &exp, &mnt, op, mode, &flags);
1990
1991    if (exp == 255 && (uint32_t)(mnt << 9)) {
1992        if (fpscr.dn) {
1993            result = fp64_defaultNaN();
1994        } else {
1995            result = fp64_FPConvertNaN_32(op);
1996        }
1997        if (!(mnt >> 22 & 1)) {
1998            flags |= FPLIB_IOC;
1999        }
2000    } else if (exp == 255) {
2001        result = fp64_infinity(sgn);
2002    } else if (!mnt) {
2003        result = fp64_zero(sgn);
2004    } else {
2005        mnt = fp32_normalise(mnt, &exp);
2006        result = fp64_pack(sgn, exp - 127 + 1023 + 8, (uint64_t)mnt << 21);
2007    }
2008
2009    set_fpscr0(fpscr, flags);
2010
2011    return result;
2012}
2013
2014template <>
2015uint32_t
2016fplibMulAdd(uint32_t addend, uint32_t op1, uint32_t op2, FPSCR &fpscr)
2017{
2018    int flags = 0;
2019    uint32_t result = fp32_muladd(addend, op1, op2, 0, modeConv(fpscr), &flags);
2020    set_fpscr0(fpscr, flags);
2021    return result;
2022}
2023
2024template <>
2025uint64_t
2026fplibMulAdd(uint64_t addend, uint64_t op1, uint64_t op2, FPSCR &fpscr)
2027{
2028    int flags = 0;
2029    uint64_t result = fp64_muladd(addend, op1, op2, 0, modeConv(fpscr), &flags);
2030    set_fpscr0(fpscr, flags);
2031    return result;
2032}
2033
2034template <>
2035uint32_t
2036fplibDiv(uint32_t op1, uint32_t op2, FPSCR &fpscr)
2037{
2038    int flags = 0;
2039    uint32_t result = fp32_div(op1, op2, modeConv(fpscr), &flags);
2040    set_fpscr0(fpscr, flags);
2041    return result;
2042}
2043
2044template <>
2045uint64_t
2046fplibDiv(uint64_t op1, uint64_t op2, FPSCR &fpscr)
2047{
2048    int flags = 0;
2049    uint64_t result = fp64_div(op1, op2, modeConv(fpscr), &flags);
2050    set_fpscr0(fpscr, flags);
2051    return result;
2052}
2053
2054static uint32_t
2055fp32_repack(int sgn, int exp, uint32_t mnt)
2056{
2057    return fp32_pack(sgn, mnt >> 23 ? exp : 0, mnt);
2058}
2059
2060static uint64_t
2061fp64_repack(int sgn, int exp, uint64_t mnt)
2062{
2063    return fp64_pack(sgn, mnt >> 52 ? exp : 0, mnt);
2064}
2065
2066static void
2067fp32_minmaxnum(uint32_t *op1, uint32_t *op2, int sgn)
2068{
2069    // Treat a single quiet-NaN as +Infinity/-Infinity
2070    if (!((uint32_t)~(*op1 << 1) >> 23) && (uint32_t)~(*op2 << 1) >> 23)
2071        *op1 = fp32_infinity(sgn);
2072    if (!((uint32_t)~(*op2 << 1) >> 23) && (uint32_t)~(*op1 << 1) >> 23)
2073        *op2 = fp32_infinity(sgn);
2074}
2075
2076static void
2077fp64_minmaxnum(uint64_t *op1, uint64_t *op2, int sgn)
2078{
2079    // Treat a single quiet-NaN as +Infinity/-Infinity
2080    if (!((uint64_t)~(*op1 << 1) >> 52) && (uint64_t)~(*op2 << 1) >> 52)
2081        *op1 = fp64_infinity(sgn);
2082    if (!((uint64_t)~(*op2 << 1) >> 52) && (uint64_t)~(*op1 << 1) >> 52)
2083        *op2 = fp64_infinity(sgn);
2084}
2085
2086template <>
2087uint32_t
2088fplibMax(uint32_t op1, uint32_t op2, FPSCR &fpscr)
2089{
2090    int mode = modeConv(fpscr);
2091    int flags = 0;
2092    int sgn1, exp1, sgn2, exp2;
2093    uint32_t mnt1, mnt2, x, result;
2094
2095    fp32_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
2096    fp32_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
2097
2098    if ((x = fp32_process_NaNs(op1, op2, mode, &flags))) {
2099        result = x;
2100    } else {
2101        result = ((sgn1 != sgn2 ? sgn2 : sgn1 ^ (op1 > op2)) ?
2102                  fp32_repack(sgn1, exp1, mnt1) :
2103                  fp32_repack(sgn2, exp2, mnt2));
2104    }
2105    set_fpscr0(fpscr, flags);
2106    return result;
2107}
2108
2109template <>
2110uint64_t
2111fplibMax(uint64_t op1, uint64_t op2, FPSCR &fpscr)
2112{
2113    int mode = modeConv(fpscr);
2114    int flags = 0;
2115    int sgn1, exp1, sgn2, exp2;
2116    uint64_t mnt1, mnt2, x, result;
2117
2118    fp64_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
2119    fp64_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
2120
2121    if ((x = fp64_process_NaNs(op1, op2, mode, &flags))) {
2122        result = x;
2123    } else {
2124        result = ((sgn1 != sgn2 ? sgn2 : sgn1 ^ (op1 > op2)) ?
2125                  fp64_repack(sgn1, exp1, mnt1) :
2126                  fp64_repack(sgn2, exp2, mnt2));
2127    }
2128    set_fpscr0(fpscr, flags);
2129    return result;
2130}
2131
2132template <>
2133uint32_t
2134fplibMaxNum(uint32_t op1, uint32_t op2, FPSCR &fpscr)
2135{
2136    fp32_minmaxnum(&op1, &op2, 1);
2137    return fplibMax<uint32_t>(op1, op2, fpscr);
2138}
2139
2140template <>
2141uint64_t
2142fplibMaxNum(uint64_t op1, uint64_t op2, FPSCR &fpscr)
2143{
2144    fp64_minmaxnum(&op1, &op2, 1);
2145    return fplibMax<uint64_t>(op1, op2, fpscr);
2146}
2147
2148template <>
2149uint32_t
2150fplibMin(uint32_t op1, uint32_t op2, FPSCR &fpscr)
2151{
2152    int mode = modeConv(fpscr);
2153    int flags = 0;
2154    int sgn1, exp1, sgn2, exp2;
2155    uint32_t mnt1, mnt2, x, result;
2156
2157    fp32_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
2158    fp32_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
2159
2160    if ((x = fp32_process_NaNs(op1, op2, mode, &flags))) {
2161        result = x;
2162    } else {
2163        result = ((sgn1 != sgn2 ? sgn1 : sgn1 ^ (op1 < op2)) ?
2164                  fp32_repack(sgn1, exp1, mnt1) :
2165                  fp32_repack(sgn2, exp2, mnt2));
2166    }
2167    set_fpscr0(fpscr, flags);
2168    return result;
2169}
2170
2171template <>
2172uint64_t
2173fplibMin(uint64_t op1, uint64_t op2, FPSCR &fpscr)
2174{
2175    int mode = modeConv(fpscr);
2176    int flags = 0;
2177    int sgn1, exp1, sgn2, exp2;
2178    uint64_t mnt1, mnt2, x, result;
2179
2180    fp64_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
2181    fp64_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
2182
2183    if ((x = fp64_process_NaNs(op1, op2, mode, &flags))) {
2184        result = x;
2185    } else {
2186        result = ((sgn1 != sgn2 ? sgn1 : sgn1 ^ (op1 < op2)) ?
2187                  fp64_repack(sgn1, exp1, mnt1) :
2188                  fp64_repack(sgn2, exp2, mnt2));
2189    }
2190    set_fpscr0(fpscr, flags);
2191    return result;
2192}
2193
2194template <>
2195uint32_t
2196fplibMinNum(uint32_t op1, uint32_t op2, FPSCR &fpscr)
2197{
2198    fp32_minmaxnum(&op1, &op2, 0);
2199    return fplibMin<uint32_t>(op1, op2, fpscr);
2200}
2201
2202template <>
2203uint64_t
2204fplibMinNum(uint64_t op1, uint64_t op2, FPSCR &fpscr)
2205{
2206    fp64_minmaxnum(&op1, &op2, 0);
2207    return fplibMin<uint64_t>(op1, op2, fpscr);
2208}
2209
2210template <>
2211uint32_t
2212fplibMul(uint32_t op1, uint32_t op2, FPSCR &fpscr)
2213{
2214    int flags = 0;
2215    uint32_t result = fp32_mul(op1, op2, modeConv(fpscr), &flags);
2216    set_fpscr0(fpscr, flags);
2217    return result;
2218}
2219
2220template <>
2221uint64_t
2222fplibMul(uint64_t op1, uint64_t op2, FPSCR &fpscr)
2223{
2224    int flags = 0;
2225    uint64_t result = fp64_mul(op1, op2, modeConv(fpscr), &flags);
2226    set_fpscr0(fpscr, flags);
2227    return result;
2228}
2229
2230template <>
2231uint32_t
2232fplibMulX(uint32_t op1, uint32_t op2, FPSCR &fpscr)
2233{
2234    int mode = modeConv(fpscr);
2235    int flags = 0;
2236    int sgn1, exp1, sgn2, exp2;
2237    uint32_t mnt1, mnt2, result;
2238
2239    fp32_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
2240    fp32_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
2241
2242    result = fp32_process_NaNs(op1, op2, mode, &flags);
2243    if (!result) {
2244        if ((exp1 == 255 && !mnt2) || (exp2 == 255 && !mnt1)) {
2245            result = fp32_FPTwo(sgn1 ^ sgn2);
2246        } else if (exp1 == 255 || exp2 == 255) {
2247            result = fp32_infinity(sgn1 ^ sgn2);
2248        } else if (!mnt1 || !mnt2) {
2249            result = fp32_zero(sgn1 ^ sgn2);
2250        } else {
2251            result = fp32_mul(op1, op2, mode, &flags);
2252        }
2253    }
2254
2255    set_fpscr0(fpscr, flags);
2256
2257    return result;
2258}
2259
2260template <>
2261uint64_t
2262fplibMulX(uint64_t op1, uint64_t op2, FPSCR &fpscr)
2263{
2264    int mode = modeConv(fpscr);
2265    int flags = 0;
2266    int sgn1, exp1, sgn2, exp2;
2267    uint64_t mnt1, mnt2, result;
2268
2269    fp64_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
2270    fp64_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
2271
2272    result = fp64_process_NaNs(op1, op2, mode, &flags);
2273    if (!result) {
2274        if ((exp1 == 2047 && !mnt2) || (exp2 == 2047 && !mnt1)) {
2275            result = fp64_FPTwo(sgn1 ^ sgn2);
2276        } else if (exp1 == 2047 || exp2 == 2047) {
2277            result = fp64_infinity(sgn1 ^ sgn2);
2278        } else if (!mnt1 || !mnt2) {
2279            result = fp64_zero(sgn1 ^ sgn2);
2280        } else {
2281            result = fp64_mul(op1, op2, mode, &flags);
2282        }
2283    }
2284
2285    set_fpscr0(fpscr, flags);
2286
2287    return result;
2288}
2289
2290template <>
2291uint32_t
2292fplibNeg(uint32_t op)
2293{
2294    return op ^ (uint32_t)1 << 31;
2295}
2296
2297template <>
2298uint64_t
2299fplibNeg(uint64_t op)
2300{
2301    return op ^ (uint64_t)1 << 63;
2302}
2303
2304static const uint8_t recip_sqrt_estimate[256] = {
2305    255, 253, 251, 249, 247, 245, 243, 242, 240, 238, 236, 234, 233, 231, 229, 228,
2306    226, 224, 223, 221, 219, 218, 216, 215, 213, 212, 210, 209, 207, 206, 204, 203,
2307    201, 200, 198, 197, 196, 194, 193, 192, 190, 189, 188, 186, 185, 184, 183, 181,
2308    180, 179, 178, 176, 175, 174, 173, 172, 170, 169, 168, 167, 166, 165, 164, 163,
2309    162, 160, 159, 158, 157, 156, 155, 154, 153, 152, 151, 150, 149, 148, 147, 146,
2310    145, 144, 143, 142, 141, 140, 140, 139, 138, 137, 136, 135, 134, 133, 132, 131,
2311    131, 130, 129, 128, 127, 126, 126, 125, 124, 123, 122, 121, 121, 120, 119, 118,
2312    118, 117, 116, 115, 114, 114, 113, 112, 111, 111, 110, 109, 109, 108, 107, 106,
2313    105, 104, 103, 101, 100,  99,  97,  96,  95,  93,  92,  91,  90,  88,  87,  86,
2314    85,  84,  82,  81,  80,  79,  78,  77,  76,  75,  74,  72,  71,  70,  69,  68,
2315    67,  66,  65,  64,  63,  62,  61,  60,  60,  59,  58,  57,  56,  55,  54,  53,
2316    52,  51,  51,  50,  49,  48,  47,  46,  46,  45,  44,  43,  42,  42,  41,  40,
2317    39,  38,  38,  37,  36,  35,  35,  34,  33,  33,  32,  31,  30,  30,  29,  28,
2318    28,  27,  26,  26,  25,  24,  24,  23,  22,  22,  21,  20,  20,  19,  19,  18,
2319    17,  17,  16,  16,  15,  14,  14,  13,  13,  12,  11,  11,  10,  10,   9,   9,
2320    8,   8,   7,   6,   6,   5,   5,   4,   4,   3,   3,   2,   2,   1,   1,   0
2321};
2322
2323template <>
2324uint32_t
2325fplibRSqrtEstimate(uint32_t op, FPSCR &fpscr)
2326{
2327    int mode = modeConv(fpscr);
2328    int flags = 0;
2329    int sgn, exp;
2330    uint32_t mnt, result;
2331
2332    fp32_unpack(&sgn, &exp, &mnt, op, mode, &flags);
2333
2334    if (exp == 255 && (uint32_t)(mnt << 9)) {
2335        result = fp32_process_NaN(op, mode, &flags);
2336    } else if (!mnt) {
2337        result = fp32_infinity(sgn);
2338        flags |= FPLIB_DZC;
2339    } else if (sgn) {
2340        result = fp32_defaultNaN();
2341        flags |= FPLIB_IOC;
2342    } else if (exp == 255) {
2343        result = fp32_zero(0);
2344    } else {
2345        exp += 8;
2346        mnt = fp32_normalise(mnt, &exp);
2347        mnt = recip_sqrt_estimate[(~exp & 1) << 7 | (mnt >> 24 & 127)];
2348        result = fp32_pack(0, (380 - exp) >> 1, mnt << 15);
2349    }
2350
2351    set_fpscr0(fpscr, flags);
2352
2353    return result;
2354}
2355
2356template <>
2357uint64_t
2358fplibRSqrtEstimate(uint64_t op, FPSCR &fpscr)
2359{
2360    int mode = modeConv(fpscr);
2361    int flags = 0;
2362    int sgn, exp;
2363    uint64_t mnt, result;
2364
2365    fp64_unpack(&sgn, &exp, &mnt, op, mode, &flags);
2366
2367    if (exp == 2047 && (uint64_t)(mnt << 12)) {
2368        result = fp64_process_NaN(op, mode, &flags);
2369    } else if (!mnt) {
2370        result = fp64_infinity(sgn);
2371        flags |= FPLIB_DZC;
2372    } else if (sgn) {
2373        result = fp64_defaultNaN();
2374        flags |= FPLIB_IOC;
2375    } else if (exp == 2047) {
2376        result = fp32_zero(0);
2377    } else {
2378        exp += 11;
2379        mnt = fp64_normalise(mnt, &exp);
2380        mnt = recip_sqrt_estimate[(~exp & 1) << 7 | (mnt >> 56 & 127)];
2381        result = fp64_pack(0, (3068 - exp) >> 1, mnt << 44);
2382    }
2383
2384    set_fpscr0(fpscr, flags);
2385
2386    return result;
2387}
2388
2389template <>
2390uint32_t
2391fplibRSqrtStepFused(uint32_t op1, uint32_t op2, FPSCR &fpscr)
2392{
2393    int mode = modeConv(fpscr);
2394    int flags = 0;
2395    int sgn1, exp1, sgn2, exp2;
2396    uint32_t mnt1, mnt2, result;
2397
2398    op1 = fplibNeg<uint32_t>(op1);
2399    fp32_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
2400    fp32_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
2401
2402    result = fp32_process_NaNs(op1, op2, mode, &flags);
2403    if (!result) {
2404        if ((exp1 == 255 && !mnt2) || (exp2 == 255 && !mnt1)) {
2405            result = fp32_FPOnePointFive(0);
2406        } else if (exp1 == 255 || exp2 == 255) {
2407            result = fp32_infinity(sgn1 ^ sgn2);
2408        } else {
2409            result = fp32_muladd(fp32_FPThree(0), op1, op2, -1, mode, &flags);
2410        }
2411    }
2412
2413    set_fpscr0(fpscr, flags);
2414
2415    return result;
2416}
2417
2418template <>
2419uint64_t
2420fplibRSqrtStepFused(uint64_t op1, uint64_t op2, FPSCR &fpscr)
2421{
2422    int mode = modeConv(fpscr);
2423    int flags = 0;
2424    int sgn1, exp1, sgn2, exp2;
2425    uint64_t mnt1, mnt2, result;
2426
2427    op1 = fplibNeg<uint64_t>(op1);
2428    fp64_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
2429    fp64_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
2430
2431    result = fp64_process_NaNs(op1, op2, mode, &flags);
2432    if (!result) {
2433        if ((exp1 == 2047 && !mnt2) || (exp2 == 2047 && !mnt1)) {
2434            result = fp64_FPOnePointFive(0);
2435        } else if (exp1 == 2047 || exp2 == 2047) {
2436            result = fp64_infinity(sgn1 ^ sgn2);
2437        } else {
2438            result = fp64_muladd(fp64_FPThree(0), op1, op2, -1, mode, &flags);
2439        }
2440    }
2441
2442    set_fpscr0(fpscr, flags);
2443
2444    return result;
2445}
2446
2447template <>
2448uint32_t
2449fplibRecipStepFused(uint32_t op1, uint32_t op2, FPSCR &fpscr)
2450{
2451    int mode = modeConv(fpscr);
2452    int flags = 0;
2453    int sgn1, exp1, sgn2, exp2;
2454    uint32_t mnt1, mnt2, result;
2455
2456    op1 = fplibNeg<uint32_t>(op1);
2457    fp32_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
2458    fp32_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
2459
2460    result = fp32_process_NaNs(op1, op2, mode, &flags);
2461    if (!result) {
2462        if ((exp1 == 255 && !mnt2) || (exp2 == 255 && !mnt1)) {
2463            result = fp32_FPTwo(0);
2464        } else if (exp1 == 255 || exp2 == 255) {
2465            result = fp32_infinity(sgn1 ^ sgn2);
2466        } else {
2467            result = fp32_muladd(fp32_FPTwo(0), op1, op2, 0, mode, &flags);
2468        }
2469    }
2470
2471    set_fpscr0(fpscr, flags);
2472
2473    return result;
2474}
2475
2476template <>
2477uint32_t
2478fplibRecipEstimate(uint32_t op, FPSCR &fpscr)
2479{
2480    int mode = modeConv(fpscr);
2481    int flags = 0;
2482    int sgn, exp;
2483    uint32_t mnt, result;
2484
2485    fp32_unpack(&sgn, &exp, &mnt, op, mode, &flags);
2486
2487    if (exp == 255 && (uint32_t)(mnt << 9)) {
2488        result = fp32_process_NaN(op, mode, &flags);
2489    } else if (exp == 255) {
2490        result = fp32_zero(sgn);
2491    } else if (!mnt) {
2492        result = fp32_infinity(sgn);
2493        flags |= FPLIB_DZC;
2494    } else if (!((uint32_t)(op << 1) >> 22)) {
2495        bool overflow_to_inf = false;
2496        switch (FPCRRounding(fpscr)) {
2497          case FPRounding_TIEEVEN:
2498            overflow_to_inf = true;
2499            break;
2500          case FPRounding_POSINF:
2501            overflow_to_inf = !sgn;
2502            break;
2503          case FPRounding_NEGINF:
2504            overflow_to_inf = sgn;
2505            break;
2506          case FPRounding_ZERO:
2507            overflow_to_inf = false;
2508            break;
2509          default:
2510            assert(0);
2511        }
2512        result = overflow_to_inf ? fp32_infinity(sgn) : fp32_max_normal(sgn);
2513        flags |= FPLIB_OFC | FPLIB_IXC;
2514    } else if (fpscr.fz && exp >= 253) {
2515        result = fp32_zero(sgn);
2516        flags |= FPLIB_UFC;
2517    } else {
2518        exp += 8;
2519        mnt = fp32_normalise(mnt, &exp);
2520        int result_exp = 253 - exp;
2521        uint32_t fraction = (((uint32_t)1 << 19) / (mnt >> 22 | 1) + 1) >> 1;
2522        fraction <<= 15;
2523        if (result_exp == 0) {
2524            fraction >>= 1;
2525        } else if (result_exp == -1) {
2526            fraction >>= 2;
2527            result_exp = 0;
2528        }
2529        result = fp32_pack(sgn, result_exp, fraction);
2530    }
2531
2532    set_fpscr0(fpscr, flags);
2533
2534    return result;
2535}
2536
2537template <>
2538uint64_t
2539fplibRecipEstimate(uint64_t op, FPSCR &fpscr)
2540{
2541    int mode = modeConv(fpscr);
2542    int flags = 0;
2543    int sgn, exp;
2544    uint64_t mnt, result;
2545
2546    fp64_unpack(&sgn, &exp, &mnt, op, mode, &flags);
2547
2548    if (exp == 2047 && (uint64_t)(mnt << 12)) {
2549        result = fp64_process_NaN(op, mode, &flags);
2550    } else if (exp == 2047) {
2551        result = fp64_zero(sgn);
2552    } else if (!mnt) {
2553        result = fp64_infinity(sgn);
2554        flags |= FPLIB_DZC;
2555    } else if (!((uint64_t)(op << 1) >> 51)) {
2556        bool overflow_to_inf = false;
2557        switch (FPCRRounding(fpscr)) {
2558          case FPRounding_TIEEVEN:
2559            overflow_to_inf = true;
2560            break;
2561          case FPRounding_POSINF:
2562            overflow_to_inf = !sgn;
2563            break;
2564          case FPRounding_NEGINF:
2565            overflow_to_inf = sgn;
2566            break;
2567          case FPRounding_ZERO:
2568            overflow_to_inf = false;
2569            break;
2570          default:
2571            assert(0);
2572        }
2573        result = overflow_to_inf ? fp64_infinity(sgn) : fp64_max_normal(sgn);
2574        flags |= FPLIB_OFC | FPLIB_IXC;
2575    } else if (fpscr.fz && exp >= 2045) {
2576        result = fp64_zero(sgn);
2577        flags |= FPLIB_UFC;
2578    } else {
2579        exp += 11;
2580        mnt = fp64_normalise(mnt, &exp);
2581        int result_exp = 2045 - exp;
2582        uint64_t fraction = (((uint32_t)1 << 19) / (mnt >> 54 | 1) + 1) >> 1;
2583        fraction <<= 44;
2584        if (result_exp == 0) {
2585            fraction >>= 1;
2586        } else if (result_exp == -1) {
2587            fraction >>= 2;
2588            result_exp = 0;
2589        }
2590        result = fp64_pack(sgn, result_exp, fraction);
2591    }
2592
2593    set_fpscr0(fpscr, flags);
2594
2595    return result;
2596}
2597
2598template <>
2599uint64_t
2600fplibRecipStepFused(uint64_t op1, uint64_t op2, FPSCR &fpscr)
2601{
2602    int mode = modeConv(fpscr);
2603    int flags = 0;
2604    int sgn1, exp1, sgn2, exp2;
2605    uint64_t mnt1, mnt2, result;
2606
2607    op1 = fplibNeg<uint64_t>(op1);
2608    fp64_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
2609    fp64_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
2610
2611    result = fp64_process_NaNs(op1, op2, mode, &flags);
2612    if (!result) {
2613        if ((exp1 == 2047 && !mnt2) || (exp2 == 2047 && !mnt1)) {
2614            result = fp64_FPTwo(0);
2615        } else if (exp1 == 2047 || exp2 == 2047) {
2616            result = fp64_infinity(sgn1 ^ sgn2);
2617        } else {
2618            result = fp64_muladd(fp64_FPTwo(0), op1, op2, 0, mode, &flags);
2619        }
2620    }
2621
2622    set_fpscr0(fpscr, flags);
2623
2624    return result;
2625}
2626
2627template <>
2628uint32_t
2629fplibRecpX(uint32_t op, FPSCR &fpscr)
2630{
2631    int mode = modeConv(fpscr);
2632    int flags = 0;
2633    int sgn, exp;
2634    uint32_t mnt, result;
2635
2636    fp32_unpack(&sgn, &exp, &mnt, op, mode, &flags);
2637
2638    if (exp == 255 && (uint32_t)(mnt << 9)) {
2639        result = fp32_process_NaN(op, mode, &flags);
2640    }
2641    else {
2642        if (!mnt) { // Zero and denormals
2643            result = fp32_pack(sgn, 254, 0);
2644        } else { // Infinities and normals
2645            result = fp32_pack(sgn, exp ^ 255, 0);
2646        }
2647    }
2648
2649    set_fpscr0(fpscr, flags);
2650
2651    return result;
2652}
2653
2654template <>
2655uint64_t
2656fplibRecpX(uint64_t op, FPSCR &fpscr)
2657{
2658    int mode = modeConv(fpscr);
2659    int flags = 0;
2660    int sgn, exp;
2661    uint64_t mnt, result;
2662
2663    fp64_unpack(&sgn, &exp, &mnt, op, mode, &flags);
2664
2665    if (exp == 2047 && (uint64_t)(mnt << 12)) {
2666        result = fp64_process_NaN(op, mode, &flags);
2667    }
2668    else {
2669        if (!mnt) { // Zero and denormals
2670            result = fp64_pack(sgn, 2046, 0);
2671        } else { // Infinities and normals
2672            result = fp64_pack(sgn, exp ^ 2047, 0);
2673        }
2674    }
2675
2676    set_fpscr0(fpscr, flags);
2677
2678    return result;
2679}
2680
2681template <>
2682uint32_t
2683fplibRoundInt(uint32_t op, FPRounding rounding, bool exact, FPSCR &fpscr)
2684{
2685    int mode = modeConv(fpscr);
2686    int flags = 0;
2687    int sgn, exp;
2688    uint32_t mnt, result;
2689
2690    // Unpack using FPCR to determine if subnormals are flushed-to-zero:
2691    fp32_unpack(&sgn, &exp, &mnt, op, mode, &flags);
2692
2693    // Handle NaNs, infinities and zeroes:
2694    if (exp == 255 && (uint32_t)(mnt << 9)) {
2695        result = fp32_process_NaN(op, mode, &flags);
2696    } else if (exp == 255) {
2697        result = fp32_infinity(sgn);
2698    } else if (!mnt) {
2699        result = fp32_zero(sgn);
2700    } else if (exp >= 150) {
2701        // There are no fractional bits
2702        result = op;
2703    } else {
2704        // Truncate towards zero:
2705        uint32_t x = 150 - exp >= 32 ? 0 : mnt >> (150 - exp);
2706        int err = exp < 118 ? 1 :
2707            (mnt << 1 >> (149 - exp) & 3) | (mnt << 2 << (exp - 118) != 0);
2708        switch (rounding) {
2709          case FPRounding_TIEEVEN:
2710            x += (err == 3 || (err == 2 && (x & 1)));
2711            break;
2712          case FPRounding_POSINF:
2713            x += err && !sgn;
2714            break;
2715          case FPRounding_NEGINF:
2716            x += err && sgn;
2717            break;
2718          case FPRounding_ZERO:
2719            break;
2720          case FPRounding_TIEAWAY:
2721            x += err >> 1;
2722            break;
2723          default:
2724            assert(0);
2725        }
2726
2727        if (x == 0) {
2728            result = fp32_zero(sgn);
2729        } else {
2730            exp = 150;
2731            mnt = fp32_normalise(x, &exp);
2732            result = fp32_pack(sgn, exp + 8, mnt >> 8);
2733        }
2734
2735        if (err && exact)
2736            flags |= FPLIB_IXC;
2737    }
2738
2739    set_fpscr0(fpscr, flags);
2740
2741    return result;
2742}
2743
2744template <>
2745uint64_t
2746fplibRoundInt(uint64_t op, FPRounding rounding, bool exact, FPSCR &fpscr)
2747{
2748    int mode = modeConv(fpscr);
2749    int flags = 0;
2750    int sgn, exp;
2751    uint64_t mnt, result;
2752
2753    // Unpack using FPCR to determine if subnormals are flushed-to-zero:
2754    fp64_unpack(&sgn, &exp, &mnt, op, mode, &flags);
2755
2756    // Handle NaNs, infinities and zeroes:
2757    if (exp == 2047 && (uint64_t)(mnt << 12)) {
2758        result = fp64_process_NaN(op, mode, &flags);
2759    } else if (exp == 2047) {
2760        result = fp64_infinity(sgn);
2761    } else if (!mnt) {
2762        result = fp64_zero(sgn);
2763    } else if (exp >= 1075) {
2764        // There are no fractional bits
2765        result = op;
2766    } else {
2767        // Truncate towards zero:
2768        uint64_t x = 1075 - exp >= 64 ? 0 : mnt >> (1075 - exp);
2769        int err = exp < 1011 ? 1 :
2770            (mnt << 1 >> (1074 - exp) & 3) | (mnt << 2 << (exp - 1011) != 0);
2771        switch (rounding) {
2772          case FPRounding_TIEEVEN:
2773            x += (err == 3 || (err == 2 && (x & 1)));
2774            break;
2775          case FPRounding_POSINF:
2776            x += err && !sgn;
2777            break;
2778          case FPRounding_NEGINF:
2779            x += err && sgn;
2780            break;
2781          case FPRounding_ZERO:
2782            break;
2783          case FPRounding_TIEAWAY:
2784            x += err >> 1;
2785            break;
2786          default:
2787            assert(0);
2788        }
2789
2790        if (x == 0) {
2791            result = fp64_zero(sgn);
2792        } else {
2793            exp = 1075;
2794            mnt = fp64_normalise(x, &exp);
2795            result = fp64_pack(sgn, exp + 11, mnt >> 11);
2796        }
2797
2798        if (err && exact)
2799            flags |= FPLIB_IXC;
2800    }
2801
2802    set_fpscr0(fpscr, flags);
2803
2804    return result;
2805}
2806
2807template <>
2808uint32_t
2809fplibSqrt(uint32_t op, FPSCR &fpscr)
2810{
2811    int flags = 0;
2812    uint32_t result = fp32_sqrt(op, modeConv(fpscr), &flags);
2813    set_fpscr0(fpscr, flags);
2814    return result;
2815}
2816
2817template <>
2818uint64_t
2819fplibSqrt(uint64_t op, FPSCR &fpscr)
2820{
2821    int flags = 0;
2822    uint64_t result = fp64_sqrt(op, modeConv(fpscr), &flags);
2823    set_fpscr0(fpscr, flags);
2824    return result;
2825}
2826
2827template <>
2828uint32_t
2829fplibSub(uint32_t op1, uint32_t op2, FPSCR &fpscr)
2830{
2831    int flags = 0;
2832    uint32_t result = fp32_add(op1, op2, 1, modeConv(fpscr), &flags);
2833    set_fpscr0(fpscr, flags);
2834    return result;
2835}
2836
2837template <>
2838uint64_t
2839fplibSub(uint64_t op1, uint64_t op2, FPSCR &fpscr)
2840{
2841    int flags = 0;
2842    uint64_t result = fp64_add(op1, op2, 1, modeConv(fpscr), &flags);
2843    set_fpscr0(fpscr, flags);
2844    return result;
2845}
2846
2847static uint64_t
2848FPToFixed_64(int sgn, int exp, uint64_t mnt, bool u, FPRounding rounding,
2849             int *flags)
2850{
2851    uint64_t x;
2852    int err;
2853
2854    if (exp > 1023 + 63) {
2855        *flags = FPLIB_IOC;
2856        return ((uint64_t)!u << 63) - !sgn;
2857    }
2858
2859    x = lsr64(mnt << 11, 1023 + 63 - exp);
2860    err = (exp > 1023 + 63 - 2 ? 0 :
2861           (lsr64(mnt << 11, 1023 + 63 - 2 - exp) & 3) |
2862           !!(mnt << 11 & (lsl64(1, 1023 + 63 - 2 - exp) - 1)));
2863
2864    switch (rounding) {
2865      case FPRounding_TIEEVEN:
2866        x += (err == 3 || (err == 2 && (x & 1)));
2867        break;
2868      case FPRounding_POSINF:
2869        x += err && !sgn;
2870        break;
2871      case FPRounding_NEGINF:
2872        x += err && sgn;
2873        break;
2874      case FPRounding_ZERO:
2875        break;
2876      case FPRounding_TIEAWAY:
2877        x += err >> 1;
2878        break;
2879      default:
2880        assert(0);
2881    }
2882
2883    if (u ? sgn && x : x > ((uint64_t)1 << 63) - !sgn) {
2884        *flags = FPLIB_IOC;
2885        return ((uint64_t)!u << 63) - !sgn;
2886    }
2887
2888    if (err) {
2889        *flags = FPLIB_IXC;
2890    }
2891
2892    return sgn ? -x : x;
2893}
2894
2895static uint32_t
2896FPToFixed_32(int sgn, int exp, uint64_t mnt, bool u, FPRounding rounding,
2897             int *flags)
2898{
2899    uint64_t x = FPToFixed_64(sgn, exp, mnt, u, rounding, flags);
2900    if (u ? x >= (uint64_t)1 << 32 :
2901        !(x < (uint64_t)1 << 31 ||
2902          (uint64_t)-x <= (uint64_t)1 << 31)) {
2903        *flags = FPLIB_IOC;
2904        x = ((uint32_t)!u << 31) - !sgn;
2905    }
2906    return x;
2907}
2908
2909template <>
2910uint32_t
2911fplibFPToFixed(uint32_t op, int fbits, bool u, FPRounding rounding, FPSCR &fpscr)
2912{
2913    int flags = 0;
2914    int sgn, exp;
2915    uint32_t mnt, result;
2916
2917    // Unpack using FPCR to determine if subnormals are flushed-to-zero:
2918    fp32_unpack(&sgn, &exp, &mnt, op, modeConv(fpscr), &flags);
2919
2920    // If NaN, set cumulative flag or take exception:
2921    if (exp == 255 && (uint32_t)(mnt << 9)) {
2922        flags = FPLIB_IOC;
2923        result = 0;
2924    } else {
2925        result = FPToFixed_32(sgn, exp + 1023 - 127 + fbits,
2926                              (uint64_t)mnt << (52 - 23), u, rounding, &flags);
2927    }
2928
2929    set_fpscr0(fpscr, flags);
2930
2931    return result;
2932}
2933
2934template <>
2935uint32_t
2936fplibFPToFixed(uint64_t op, int fbits, bool u, FPRounding rounding, FPSCR &fpscr)
2937{
2938    int flags = 0;
2939    int sgn, exp;
2940    uint64_t mnt;
2941    uint32_t result;
2942
2943    // Unpack using FPCR to determine if subnormals are flushed-to-zero:
2944    fp64_unpack(&sgn, &exp, &mnt, op, modeConv(fpscr), &flags);
2945
2946    // If NaN, set cumulative flag or take exception:
2947    if (exp == 2047 && (uint64_t)(mnt << 12)) {
2948        flags = FPLIB_IOC;
2949        result = 0;
2950    } else {
2951        result = FPToFixed_32(sgn, exp + fbits, mnt, u, rounding, &flags);
2952    }
2953
2954    set_fpscr0(fpscr, flags);
2955
2956    return result;
2957}
2958
2959template <>
2960uint64_t
2961fplibFPToFixed(uint32_t op, int fbits, bool u, FPRounding rounding, FPSCR &fpscr)
2962{
2963    int flags = 0;
2964    int sgn, exp;
2965    uint32_t mnt;
2966    uint64_t result;
2967
2968    // Unpack using FPCR to determine if subnormals are flushed-to-zero:
2969    fp32_unpack(&sgn, &exp, &mnt, op, modeConv(fpscr), &flags);
2970
2971    // If NaN, set cumulative flag or take exception:
2972    if (exp == 255 && (uint32_t)(mnt << 9)) {
2973        flags = FPLIB_IOC;
2974        result = 0;
2975    } else {
2976        result = FPToFixed_64(sgn, exp + 1023 - 127 + fbits,
2977                              (uint64_t)mnt << (52 - 23), u, rounding, &flags);
2978    }
2979
2980    set_fpscr0(fpscr, flags);
2981
2982    return result;
2983}
2984
2985template <>
2986uint64_t
2987fplibFPToFixed(uint64_t op, int fbits, bool u, FPRounding rounding, FPSCR &fpscr)
2988{
2989    int flags = 0;
2990    int sgn, exp;
2991    uint64_t mnt, result;
2992
2993    // Unpack using FPCR to determine if subnormals are flushed-to-zero:
2994    fp64_unpack(&sgn, &exp, &mnt, op, modeConv(fpscr), &flags);
2995
2996    // If NaN, set cumulative flag or take exception:
2997    if (exp == 2047 && (uint64_t)(mnt << 12)) {
2998        flags = FPLIB_IOC;
2999        result = 0;
3000    } else {
3001        result = FPToFixed_64(sgn, exp + fbits, mnt, u, rounding, &flags);
3002    }
3003
3004    set_fpscr0(fpscr, flags);
3005
3006    return result;
3007}
3008
3009static uint32_t
3010fp32_cvtf(uint64_t a, int fbits, int u, int mode, int *flags)
3011{
3012    int x_sgn = !u && a >> 63;
3013    int x_exp = 190 - fbits;
3014    uint64_t x_mnt = x_sgn ? -a : a;
3015
3016    // Handle zero:
3017    if (!x_mnt) {
3018        return fp32_zero(0);
3019    }
3020
3021    // Normalise and convert to 32 bits, collapsing error into bottom bit:
3022    x_mnt = fp64_normalise(x_mnt, &x_exp);
3023    x_mnt = x_mnt >> 31 | !!(uint32_t)(x_mnt << 1);
3024
3025    return fp32_round(x_sgn, x_exp, x_mnt, mode, flags);
3026}
3027
3028static uint64_t
3029fp64_cvtf(uint64_t a, int fbits, int u, int mode, int *flags)
3030{
3031    int x_sgn = !u && a >> 63;
3032    int x_exp = 1024 + 62 - fbits;
3033    uint64_t x_mnt = x_sgn ? -a : a;
3034
3035    // Handle zero:
3036    if (!x_mnt) {
3037        return fp64_zero(0);
3038    }
3039
3040    x_mnt = fp64_normalise(x_mnt, &x_exp);
3041
3042    return fp64_round(x_sgn, x_exp, x_mnt << 1, mode, flags);
3043}
3044
3045template <>
3046uint32_t
3047fplibFixedToFP(uint64_t op, int fbits, bool u, FPRounding rounding, FPSCR &fpscr)
3048{
3049    int flags = 0;
3050    uint32_t res = fp32_cvtf(op, fbits, u,
3051                             (int)rounding | ((uint32_t)fpscr >> 22 & 12),
3052                             &flags);
3053    set_fpscr0(fpscr, flags);
3054    return res;
3055}
3056
3057template <>
3058uint64_t
3059fplibFixedToFP(uint64_t op, int fbits, bool u, FPRounding rounding, FPSCR &fpscr)
3060{
3061    int flags = 0;
3062    uint64_t res = fp64_cvtf(op, fbits, u,
3063                             (int)rounding | ((uint32_t)fpscr >> 22 & 12),
3064                             &flags);
3065    set_fpscr0(fpscr, flags);
3066    return res;
3067}
3068
3069}
3070