fplib.cc revision 13449:2f7efa89c58b 
1/*
2* Copyright (c) 2012-2013, 2017-2018 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 "base/logging.hh"
46#include "fplib.hh"
47
48namespace ArmISA
49{
50
51#define FPLIB_RN 0
52#define FPLIB_RP 1
53#define FPLIB_RM 2
54#define FPLIB_RZ 3
55#define FPLIB_FZ 4
56#define FPLIB_DN 8
57#define FPLIB_AHP 16
58#define FPLIB_FZ16 32
59
60#define FPLIB_IDC 128 // Input Denormal
61#define FPLIB_IXC 16  // Inexact
62#define FPLIB_UFC 8   // Underflow
63#define FPLIB_OFC 4   // Overflow
64#define FPLIB_DZC 2   // Division by Zero
65#define FPLIB_IOC 1   // Invalid Operation
66
67#define FP16_BITS 16
68#define FP32_BITS 32
69#define FP64_BITS 64
70
71#define FP16_EXP_BITS 5
72#define FP32_EXP_BITS 8
73#define FP64_EXP_BITS 11
74
75#define FP16_EXP_BIAS 15
76#define FP32_EXP_BIAS 127
77#define FP64_EXP_BIAS 1023
78
79#define FP16_EXP_INF ((1ULL << FP16_EXP_BITS) - 1)
80#define FP32_EXP_INF ((1ULL << FP32_EXP_BITS) - 1)
81#define FP64_EXP_INF ((1ULL << FP64_EXP_BITS) - 1)
82
83#define FP16_MANT_BITS (FP16_BITS - FP16_EXP_BITS - 1)
84#define FP32_MANT_BITS (FP32_BITS - FP32_EXP_BITS - 1)
85#define FP64_MANT_BITS (FP64_BITS - FP64_EXP_BITS - 1)
86
87#define FP16_EXP(x) ((x) >> FP16_MANT_BITS & ((1ULL << FP16_EXP_BITS) - 1))
88#define FP32_EXP(x) ((x) >> FP32_MANT_BITS & ((1ULL << FP32_EXP_BITS) - 1))
89#define FP64_EXP(x) ((x) >> FP64_MANT_BITS & ((1ULL << FP64_EXP_BITS) - 1))
90
91#define FP16_MANT(x) ((x) & ((1ULL << FP16_MANT_BITS) - 1))
92#define FP32_MANT(x) ((x) & ((1ULL << FP32_MANT_BITS) - 1))
93#define FP64_MANT(x) ((x) & ((1ULL << FP64_MANT_BITS) - 1))
94
95static inline uint16_t
96lsl16(uint16_t x, uint32_t shift)
97{
98    return shift < 16 ? x << shift : 0;
99}
100
101static inline uint16_t
102lsr16(uint16_t x, uint32_t shift)
103{
104    return shift < 16 ? x >> shift : 0;
105}
106
107static inline uint32_t
108lsl32(uint32_t x, uint32_t shift)
109{
110    return shift < 32 ? x << shift : 0;
111}
112
113static inline uint32_t
114lsr32(uint32_t x, uint32_t shift)
115{
116    return shift < 32 ? x >> shift : 0;
117}
118
119static inline uint64_t
120lsl64(uint64_t x, uint32_t shift)
121{
122    return shift < 64 ? x << shift : 0;
123}
124
125static inline uint64_t
126lsr64(uint64_t x, uint32_t shift)
127{
128    return shift < 64 ? x >> shift : 0;
129}
130
131static inline void
132lsl128(uint64_t *r0, uint64_t *r1, uint64_t x0, uint64_t x1, uint32_t shift)
133{
134    if (shift == 0) {
135        *r1 = x1;
136        *r0 = x0;
137    } else if (shift < 64) {
138        *r1 = x1 << shift | x0 >> (64 - shift);
139        *r0 = x0 << shift;
140    } else if (shift < 128) {
141        *r1 = x0 << (shift - 64);
142        *r0 = 0;
143    } else {
144        *r1 = 0;
145        *r0 = 0;
146    }
147}
148
149static inline void
150lsr128(uint64_t *r0, uint64_t *r1, uint64_t x0, uint64_t x1, uint32_t shift)
151{
152    if (shift == 0) {
153        *r1 = x1;
154        *r0 = x0;
155    } else if (shift < 64) {
156        *r0 = x0 >> shift | x1 << (64 - shift);
157        *r1 = x1 >> shift;
158    } else if (shift < 128) {
159        *r0 = x1 >> (shift - 64);
160        *r1 = 0;
161    } else {
162        *r0 = 0;
163        *r1 = 0;
164    }
165}
166
167static inline void
168mul62x62(uint64_t *x0, uint64_t *x1, uint64_t a, uint64_t b)
169{
170    uint32_t mask = ((uint32_t)1 << 31) - 1;
171    uint64_t a0 = a & mask;
172    uint64_t a1 = a >> 31 & mask;
173    uint64_t b0 = b & mask;
174    uint64_t b1 = b >> 31 & mask;
175    uint64_t p0 = a0 * b0;
176    uint64_t p2 = a1 * b1;
177    uint64_t p1 = (a0 + a1) * (b0 + b1) - p0 - p2;
178    uint64_t s0 = p0;
179    uint64_t s1 = (s0 >> 31) + p1;
180    uint64_t s2 = (s1 >> 31) + p2;
181    *x0 = (s0 & mask) | (s1 & mask) << 31 | s2 << 62;
182    *x1 = s2 >> 2;
183}
184
185static inline
186void mul64x32(uint64_t *x0, uint64_t *x1, uint64_t a, uint32_t b)
187{
188    uint64_t t0 = (uint64_t)(uint32_t)a * b;
189    uint64_t t1 = (t0 >> 32) + (a >> 32) * b;
190    *x0 = t1 << 32 | (uint32_t)t0;
191    *x1 = t1 >> 32;
192}
193
194static inline void
195add128(uint64_t *x0, uint64_t *x1, uint64_t a0, uint64_t a1, uint64_t b0,
196       uint64_t b1)
197{
198    *x0 = a0 + b0;
199    *x1 = a1 + b1 + (*x0 < a0);
200}
201
202static inline void
203sub128(uint64_t *x0, uint64_t *x1, uint64_t a0, uint64_t a1, uint64_t b0,
204       uint64_t b1)
205{
206    *x0 = a0 - b0;
207    *x1 = a1 - b1 - (*x0 > a0);
208}
209
210static inline int
211cmp128(uint64_t a0, uint64_t a1, uint64_t b0, uint64_t b1)
212{
213    return (a1 < b1 ? -1 : a1 > b1 ? 1 : a0 < b0 ? -1 : a0 > b0 ? 1 : 0);
214}
215
216static inline uint16_t
217fp16_normalise(uint16_t mnt, int *exp)
218{
219    int shift;
220
221    if (!mnt) {
222        return 0;
223    }
224
225    for (shift = 8; shift; shift >>= 1) {
226        if (!(mnt >> (16 - shift))) {
227            mnt <<= shift;
228            *exp -= shift;
229        }
230    }
231    return mnt;
232}
233
234static inline uint32_t
235fp32_normalise(uint32_t mnt, int *exp)
236{
237    int shift;
238
239    if (!mnt) {
240        return 0;
241    }
242
243    for (shift = 16; shift; shift >>= 1) {
244        if (!(mnt >> (32 - shift))) {
245            mnt <<= shift;
246            *exp -= shift;
247        }
248    }
249    return mnt;
250}
251
252static inline uint64_t
253fp64_normalise(uint64_t mnt, int *exp)
254{
255    int shift;
256
257    if (!mnt) {
258        return 0;
259    }
260
261    for (shift = 32; shift; shift >>= 1) {
262        if (!(mnt >> (64 - shift))) {
263            mnt <<= shift;
264            *exp -= shift;
265        }
266    }
267    return mnt;
268}
269
270static inline void
271fp128_normalise(uint64_t *mnt0, uint64_t *mnt1, int *exp)
272{
273    uint64_t x0 = *mnt0;
274    uint64_t x1 = *mnt1;
275    int shift;
276
277    if (!x0 && !x1) {
278        return;
279    }
280
281    if (!x1) {
282        x1 = x0;
283        x0 = 0;
284        *exp -= 64;
285    }
286
287    for (shift = 32; shift; shift >>= 1) {
288        if (!(x1 >> (64 - shift))) {
289            x1 = x1 << shift | x0 >> (64 - shift);
290            x0 <<= shift;
291            *exp -= shift;
292        }
293    }
294
295    *mnt0 = x0;
296    *mnt1 = x1;
297}
298
299static inline uint16_t
300fp16_pack(uint16_t sgn, uint16_t exp, uint16_t mnt)
301{
302    return sgn << (FP16_BITS - 1) | exp << FP16_MANT_BITS | FP16_MANT(mnt);
303}
304
305static inline uint32_t
306fp32_pack(uint32_t sgn, uint32_t exp, uint32_t mnt)
307{
308    return sgn << (FP32_BITS - 1) | exp << FP32_MANT_BITS | FP32_MANT(mnt);
309}
310
311static inline uint64_t
312fp64_pack(uint64_t sgn, uint64_t exp, uint64_t mnt)
313{
314    return sgn << (FP64_BITS - 1) | exp << FP64_MANT_BITS | FP64_MANT(mnt);
315}
316
317static inline uint16_t
318fp16_zero(int sgn)
319{
320    return fp16_pack(sgn, 0, 0);
321}
322
323static inline uint32_t
324fp32_zero(int sgn)
325{
326    return fp32_pack(sgn, 0, 0);
327}
328
329static inline uint64_t
330fp64_zero(int sgn)
331{
332    return fp64_pack(sgn, 0, 0);
333}
334
335static inline uint16_t
336fp16_max_normal(int sgn)
337{
338    return fp16_pack(sgn, FP16_EXP_INF - 1, -1);
339}
340
341static inline uint32_t
342fp32_max_normal(int sgn)
343{
344    return fp32_pack(sgn, FP32_EXP_INF - 1, -1);
345}
346
347static inline uint64_t
348fp64_max_normal(int sgn)
349{
350    return fp64_pack(sgn, FP64_EXP_INF - 1, -1);
351}
352
353static inline uint16_t
354fp16_infinity(int sgn)
355{
356    return fp16_pack(sgn, FP16_EXP_INF, 0);
357}
358
359static inline uint32_t
360fp32_infinity(int sgn)
361{
362    return fp32_pack(sgn, FP32_EXP_INF, 0);
363}
364
365static inline uint64_t
366fp64_infinity(int sgn)
367{
368    return fp64_pack(sgn, FP64_EXP_INF, 0);
369}
370
371static inline uint16_t
372fp16_defaultNaN()
373{
374    return fp16_pack(0, FP16_EXP_INF, 1ULL << (FP16_MANT_BITS - 1));
375}
376
377static inline uint32_t
378fp32_defaultNaN()
379{
380    return fp32_pack(0, FP32_EXP_INF, 1ULL << (FP32_MANT_BITS - 1));
381}
382
383static inline uint64_t
384fp64_defaultNaN()
385{
386    return fp64_pack(0, FP64_EXP_INF, 1ULL << (FP64_MANT_BITS - 1));
387}
388
389static inline void
390fp16_unpack(int *sgn, int *exp, uint16_t *mnt, uint16_t x, int mode,
391            int *flags)
392{
393    *sgn = x >> (FP16_BITS - 1);
394    *exp = FP16_EXP(x);
395    *mnt = FP16_MANT(x);
396
397    // Handle subnormals:
398    if (*exp) {
399        *mnt |= 1ULL << FP16_MANT_BITS;
400    } else {
401        ++*exp;
402        // IDC (Input Denormal) is not set in this case.
403        if (mode & FPLIB_FZ16)
404            *mnt = 0;
405    }
406}
407
408static inline void
409fp32_unpack(int *sgn, int *exp, uint32_t *mnt, uint32_t x, int mode,
410            int *flags)
411{
412    *sgn = x >> (FP32_BITS - 1);
413    *exp = FP32_EXP(x);
414    *mnt = FP32_MANT(x);
415
416    // Handle subnormals:
417    if (*exp) {
418        *mnt |= 1ULL << FP32_MANT_BITS;
419    } else {
420        ++*exp;
421        if ((mode & FPLIB_FZ) && *mnt) {
422            *flags |= FPLIB_IDC;
423            *mnt = 0;
424        }
425    }
426}
427
428static inline void
429fp64_unpack(int *sgn, int *exp, uint64_t *mnt, uint64_t x, int mode,
430            int *flags)
431{
432    *sgn = x >> (FP64_BITS - 1);
433    *exp = FP64_EXP(x);
434    *mnt = FP64_MANT(x);
435
436    // Handle subnormals:
437    if (*exp) {
438        *mnt |= 1ULL << FP64_MANT_BITS;
439    } else {
440        ++*exp;
441        if ((mode & FPLIB_FZ) && *mnt) {
442            *flags |= FPLIB_IDC;
443            *mnt = 0;
444        }
445    }
446}
447
448static inline int
449fp16_is_NaN(int exp, uint16_t mnt)
450{
451    return exp == FP16_EXP_INF && FP16_MANT(mnt);
452}
453
454static inline int
455fp32_is_NaN(int exp, uint32_t mnt)
456{
457    return exp == FP32_EXP_INF && FP32_MANT(mnt);
458}
459
460static inline int
461fp64_is_NaN(int exp, uint64_t mnt)
462{
463    return exp == FP64_EXP_INF && FP64_MANT(mnt);
464}
465
466static inline int
467fp16_is_signalling_NaN(int exp, uint16_t mnt)
468{
469    return fp16_is_NaN(exp, mnt) && !(mnt >> (FP16_MANT_BITS - 1) & 1);
470}
471
472static inline int
473fp32_is_signalling_NaN(int exp, uint32_t mnt)
474{
475    return fp32_is_NaN(exp, mnt) && !(mnt >> (FP32_MANT_BITS - 1) & 1);
476}
477
478static inline int
479fp64_is_signalling_NaN(int exp, uint64_t mnt)
480{
481    return fp64_is_NaN(exp, mnt) && !(mnt >> (FP64_MANT_BITS - 1) & 1);
482}
483
484static inline int
485fp16_is_quiet_NaN(int exp, uint16_t mnt)
486{
487    return exp == FP16_EXP_INF && (mnt >> (FP16_MANT_BITS - 1) & 1);
488}
489
490static inline int
491fp32_is_quiet_NaN(int exp, uint32_t mnt)
492{
493    return exp == FP32_EXP_INF && (mnt >> (FP32_MANT_BITS - 1) & 1);
494}
495
496static inline int
497fp64_is_quiet_NaN(int exp, uint64_t mnt)
498{
499    return exp == FP64_EXP_INF && (mnt >> (FP64_MANT_BITS - 1) & 1);
500}
501
502static inline int
503fp16_is_infinity(int exp, uint16_t mnt)
504{
505    return exp == FP16_EXP_INF && !FP16_MANT(mnt);
506}
507
508static inline int
509fp32_is_infinity(int exp, uint32_t mnt)
510{
511    return exp == FP32_EXP_INF && !FP32_MANT(mnt);
512}
513
514static inline int
515fp64_is_infinity(int exp, uint64_t mnt)
516{
517    return exp == FP64_EXP_INF && !FP64_MANT(mnt);
518}
519
520static inline uint16_t
521fp16_process_NaN(uint16_t a, int mode, int *flags)
522{
523    if (!(a >> (FP16_MANT_BITS - 1) & 1)) {
524        *flags |= FPLIB_IOC;
525        a |= 1ULL << (FP16_MANT_BITS - 1);
526    }
527    return mode & FPLIB_DN ? fp16_defaultNaN() : a;
528}
529
530static inline uint32_t
531fp32_process_NaN(uint32_t a, int mode, int *flags)
532{
533    if (!(a >> (FP32_MANT_BITS - 1) & 1)) {
534        *flags |= FPLIB_IOC;
535        a |= 1ULL << (FP32_MANT_BITS - 1);
536    }
537    return mode & FPLIB_DN ? fp32_defaultNaN() : a;
538}
539
540static inline uint64_t
541fp64_process_NaN(uint64_t a, int mode, int *flags)
542{
543    if (!(a >> (FP64_MANT_BITS - 1) & 1)) {
544        *flags |= FPLIB_IOC;
545        a |= 1ULL << (FP64_MANT_BITS - 1);
546    }
547    return mode & FPLIB_DN ? fp64_defaultNaN() : a;
548}
549
550static uint16_t
551fp16_process_NaNs(uint16_t a, uint16_t b, int mode, int *flags)
552{
553    int a_exp = FP16_EXP(a);
554    uint16_t a_mnt = FP16_MANT(a);
555    int b_exp = FP16_EXP(b);
556    uint16_t b_mnt = FP16_MANT(b);
557
558    // Handle signalling NaNs:
559    if (fp16_is_signalling_NaN(a_exp, a_mnt))
560        return fp16_process_NaN(a, mode, flags);
561    if (fp16_is_signalling_NaN(b_exp, b_mnt))
562        return fp16_process_NaN(b, mode, flags);
563
564    // Handle quiet NaNs:
565    if (fp16_is_NaN(a_exp, a_mnt))
566        return fp16_process_NaN(a, mode, flags);
567    if (fp16_is_NaN(b_exp, b_mnt))
568        return fp16_process_NaN(b, mode, flags);
569
570    return 0;
571}
572
573static uint32_t
574fp32_process_NaNs(uint32_t a, uint32_t b, int mode, int *flags)
575{
576    int a_exp = FP32_EXP(a);
577    uint32_t a_mnt = FP32_MANT(a);
578    int b_exp = FP32_EXP(b);
579    uint32_t b_mnt = FP32_MANT(b);
580
581    // Handle signalling NaNs:
582    if (fp32_is_signalling_NaN(a_exp, a_mnt))
583        return fp32_process_NaN(a, mode, flags);
584    if (fp32_is_signalling_NaN(b_exp, b_mnt))
585        return fp32_process_NaN(b, mode, flags);
586
587    // Handle quiet NaNs:
588    if (fp32_is_NaN(a_exp, a_mnt))
589        return fp32_process_NaN(a, mode, flags);
590    if (fp32_is_NaN(b_exp, b_mnt))
591        return fp32_process_NaN(b, mode, flags);
592
593    return 0;
594}
595
596static uint64_t
597fp64_process_NaNs(uint64_t a, uint64_t b, int mode, int *flags)
598{
599    int a_exp = FP64_EXP(a);
600    uint64_t a_mnt = FP64_MANT(a);
601    int b_exp = FP64_EXP(b);
602    uint64_t b_mnt = FP64_MANT(b);
603
604    // Handle signalling NaNs:
605    if (fp64_is_signalling_NaN(a_exp, a_mnt))
606        return fp64_process_NaN(a, mode, flags);
607    if (fp64_is_signalling_NaN(b_exp, b_mnt))
608        return fp64_process_NaN(b, mode, flags);
609
610    // Handle quiet NaNs:
611    if (fp64_is_NaN(a_exp, a_mnt))
612        return fp64_process_NaN(a, mode, flags);
613    if (fp64_is_NaN(b_exp, b_mnt))
614        return fp64_process_NaN(b, mode, flags);
615
616    return 0;
617}
618
619static uint16_t
620fp16_process_NaNs3(uint16_t a, uint16_t b, uint16_t c, int mode, int *flags)
621{
622    int a_exp = FP16_EXP(a);
623    uint16_t a_mnt = FP16_MANT(a);
624    int b_exp = FP16_EXP(b);
625    uint16_t b_mnt = FP16_MANT(b);
626    int c_exp = FP16_EXP(c);
627    uint16_t c_mnt = FP16_MANT(c);
628
629    // Handle signalling NaNs:
630    if (fp16_is_signalling_NaN(a_exp, a_mnt))
631        return fp16_process_NaN(a, mode, flags);
632    if (fp16_is_signalling_NaN(b_exp, b_mnt))
633        return fp16_process_NaN(b, mode, flags);
634    if (fp16_is_signalling_NaN(c_exp, c_mnt))
635        return fp16_process_NaN(c, mode, flags);
636
637    // Handle quiet NaNs:
638    if (fp16_is_NaN(a_exp, a_mnt))
639        return fp16_process_NaN(a, mode, flags);
640    if (fp16_is_NaN(b_exp, b_mnt))
641        return fp16_process_NaN(b, mode, flags);
642    if (fp16_is_NaN(c_exp, c_mnt))
643        return fp16_process_NaN(c, mode, flags);
644
645    return 0;
646}
647
648static uint32_t
649fp32_process_NaNs3(uint32_t a, uint32_t b, uint32_t c, int mode, int *flags)
650{
651    int a_exp = FP32_EXP(a);
652    uint32_t a_mnt = FP32_MANT(a);
653    int b_exp = FP32_EXP(b);
654    uint32_t b_mnt = FP32_MANT(b);
655    int c_exp = FP32_EXP(c);
656    uint32_t c_mnt = FP32_MANT(c);
657
658    // Handle signalling NaNs:
659    if (fp32_is_signalling_NaN(a_exp, a_mnt))
660        return fp32_process_NaN(a, mode, flags);
661    if (fp32_is_signalling_NaN(b_exp, b_mnt))
662        return fp32_process_NaN(b, mode, flags);
663    if (fp32_is_signalling_NaN(c_exp, c_mnt))
664        return fp32_process_NaN(c, mode, flags);
665
666    // Handle quiet NaNs:
667    if (fp32_is_NaN(a_exp, a_mnt))
668        return fp32_process_NaN(a, mode, flags);
669    if (fp32_is_NaN(b_exp, b_mnt))
670        return fp32_process_NaN(b, mode, flags);
671    if (fp32_is_NaN(c_exp, c_mnt))
672        return fp32_process_NaN(c, mode, flags);
673
674    return 0;
675}
676
677static uint64_t
678fp64_process_NaNs3(uint64_t a, uint64_t b, uint64_t c, int mode, int *flags)
679{
680    int a_exp = FP64_EXP(a);
681    uint64_t a_mnt = FP64_MANT(a);
682    int b_exp = FP64_EXP(b);
683    uint64_t b_mnt = FP64_MANT(b);
684    int c_exp = FP64_EXP(c);
685    uint64_t c_mnt = FP64_MANT(c);
686
687    // Handle signalling NaNs:
688    if (fp64_is_signalling_NaN(a_exp, a_mnt))
689        return fp64_process_NaN(a, mode, flags);
690    if (fp64_is_signalling_NaN(b_exp, b_mnt))
691        return fp64_process_NaN(b, mode, flags);
692    if (fp64_is_signalling_NaN(c_exp, c_mnt))
693        return fp64_process_NaN(c, mode, flags);
694
695    // Handle quiet NaNs:
696    if (fp64_is_NaN(a_exp, a_mnt))
697        return fp64_process_NaN(a, mode, flags);
698    if (fp64_is_NaN(b_exp, b_mnt))
699        return fp64_process_NaN(b, mode, flags);
700    if (fp64_is_NaN(c_exp, c_mnt))
701        return fp64_process_NaN(c, mode, flags);
702
703    return 0;
704}
705
706static uint16_t
707fp16_round_(int sgn, int exp, uint16_t mnt, int rm, int mode, int *flags)
708{
709    int biased_exp; // non-negative exponent value for result
710    uint16_t int_mant; // mantissa for result, less than (2 << FP16_MANT_BITS)
711    int error; // 0, 1, 2 or 3, where 2 means int_mant is wrong by exactly 0.5
712
713    assert(rm != FPRounding_TIEAWAY);
714
715    // Flush to zero:
716    if ((mode & FPLIB_FZ16) && exp < 1) {
717        *flags |= FPLIB_UFC;
718        return fp16_zero(sgn);
719    }
720
721    // The bottom FP16_EXP_BITS bits of mnt are orred together:
722    mnt = (4ULL << FP16_MANT_BITS | mnt >> (FP16_EXP_BITS - 1) |
723           ((mnt & ((1ULL << FP16_EXP_BITS) - 1)) != 0));
724
725    if (exp > 0) {
726        biased_exp = exp;
727        int_mant = mnt >> 2;
728        error = mnt & 3;
729    } else {
730        biased_exp = 0;
731        int_mant = lsr16(mnt, 3 - exp);
732        error = (lsr16(mnt, 1 - exp) & 3) | !!(mnt & (lsl16(1, 1 - exp) - 1));
733    }
734
735    if (!biased_exp && error) { // xx should also check fpscr_val<11>
736        *flags |= FPLIB_UFC;
737    }
738
739    // Round up:
740    if ((rm == FPLIB_RN && (error == 3 ||
741                            (error == 2 && (int_mant & 1)))) ||
742        (((rm == FPLIB_RP && !sgn) || (rm == FPLIB_RM && sgn)) && error)) {
743        ++int_mant;
744        if (int_mant == 1ULL << FP16_MANT_BITS) {
745            // Rounded up from denormalized to normalized
746            biased_exp = 1;
747        }
748        if (int_mant == 2ULL << FP16_MANT_BITS) {
749            // Rounded up to next exponent
750            ++biased_exp;
751            int_mant >>= 1;
752        }
753    }
754
755    // Handle rounding to odd aka Von Neumann rounding:
756    if (error && rm == FPRounding_ODD)
757        int_mant |= 1;
758
759    // Handle overflow:
760    if (!(mode & FPLIB_AHP)) {
761        if (biased_exp >= (int)FP16_EXP_INF) {
762            *flags |= FPLIB_OFC | FPLIB_IXC;
763            if (rm == FPLIB_RN || (rm == FPLIB_RP && !sgn) ||
764                (rm == FPLIB_RM && sgn)) {
765                return fp16_infinity(sgn);
766            } else {
767                return fp16_max_normal(sgn);
768            }
769        }
770    } else {
771        if (biased_exp >= (int)FP16_EXP_INF + 1) {
772            *flags |= FPLIB_IOC;
773            return fp16_pack(sgn, FP16_EXP_INF, -1);
774        }
775    }
776
777    if (error) {
778        *flags |= FPLIB_IXC;
779    }
780
781    return fp16_pack(sgn, biased_exp, int_mant);
782}
783
784static uint16_t
785fp16_round(int sgn, int exp, uint16_t mnt, int mode, int *flags)
786{
787    return fp16_round_(sgn, exp, mnt, mode & 3, mode, flags);
788}
789
790static uint32_t
791fp32_round_(int sgn, int exp, uint32_t mnt, int rm, int mode, int *flags)
792{
793    int biased_exp; // non-negative exponent value for result
794    uint32_t int_mant; // mantissa for result, less than (2 << FP32_MANT_BITS)
795    int error; // 0, 1, 2 or 3, where 2 means int_mant is wrong by exactly 0.5
796
797    assert(rm != FPRounding_TIEAWAY);
798
799    // Flush to zero:
800    if ((mode & FPLIB_FZ) && exp < 1) {
801        *flags |= FPLIB_UFC;
802        return fp32_zero(sgn);
803    }
804
805    // The bottom FP32_EXP_BITS bits of mnt are orred together:
806    mnt = (4ULL << FP32_MANT_BITS | mnt >> (FP32_EXP_BITS - 1) |
807           ((mnt & ((1ULL << FP32_EXP_BITS) - 1)) != 0));
808
809    if (exp > 0) {
810        biased_exp = exp;
811        int_mant = mnt >> 2;
812        error = mnt & 3;
813    } else {
814        biased_exp = 0;
815        int_mant = lsr32(mnt, 3 - exp);
816        error = (lsr32(mnt, 1 - exp) & 3) | !!(mnt & (lsl32(1, 1 - exp) - 1));
817    }
818
819    if (!biased_exp && error) { // xx should also check fpscr_val<11>
820        *flags |= FPLIB_UFC;
821    }
822
823    // Round up:
824    if ((rm == FPLIB_RN && (error == 3 ||
825                            (error == 2 && (int_mant & 1)))) ||
826        (((rm == FPLIB_RP && !sgn) || (rm == FPLIB_RM && sgn)) && error)) {
827        ++int_mant;
828        if (int_mant == 1ULL << FP32_MANT_BITS) {
829            // Rounded up from denormalized to normalized
830            biased_exp = 1;
831        }
832        if (int_mant == 2ULL << FP32_MANT_BITS) {
833            // Rounded up to next exponent
834            ++biased_exp;
835            int_mant >>= 1;
836        }
837    }
838
839    // Handle rounding to odd aka Von Neumann rounding:
840    if (error && rm == FPRounding_ODD)
841        int_mant |= 1;
842
843    // Handle overflow:
844    if (biased_exp >= (int)FP32_EXP_INF) {
845        *flags |= FPLIB_OFC | FPLIB_IXC;
846        if (rm == FPLIB_RN || (rm == FPLIB_RP && !sgn) ||
847            (rm == FPLIB_RM && sgn)) {
848            return fp32_infinity(sgn);
849        } else {
850            return fp32_max_normal(sgn);
851        }
852    }
853
854    if (error) {
855        *flags |= FPLIB_IXC;
856    }
857
858    return fp32_pack(sgn, biased_exp, int_mant);
859}
860
861static uint32_t
862fp32_round(int sgn, int exp, uint32_t mnt, int mode, int *flags)
863{
864    return fp32_round_(sgn, exp, mnt, mode & 3, mode, flags);
865}
866
867static uint64_t
868fp64_round_(int sgn, int exp, uint64_t mnt, int rm, int mode, int *flags)
869{
870    int biased_exp; // non-negative exponent value for result
871    uint64_t int_mant; // mantissa for result, less than (2 << FP64_MANT_BITS)
872    int error; // 0, 1, 2 or 3, where 2 means int_mant is wrong by exactly 0.5
873
874    assert(rm != FPRounding_TIEAWAY);
875
876    // Flush to zero:
877    if ((mode & FPLIB_FZ) && exp < 1) {
878        *flags |= FPLIB_UFC;
879        return fp64_zero(sgn);
880    }
881
882    // The bottom FP64_EXP_BITS bits of mnt are orred together:
883    mnt = (4ULL << FP64_MANT_BITS | mnt >> (FP64_EXP_BITS - 1) |
884           ((mnt & ((1ULL << FP64_EXP_BITS) - 1)) != 0));
885
886    if (exp > 0) {
887        biased_exp = exp;
888        int_mant = mnt >> 2;
889        error = mnt & 3;
890    } else {
891        biased_exp = 0;
892        int_mant = lsr64(mnt, 3 - exp);
893        error = (lsr64(mnt, 1 - exp) & 3) | !!(mnt & (lsl64(1, 1 - exp) - 1));
894    }
895
896    if (!biased_exp && error) { // xx should also check fpscr_val<11>
897        *flags |= FPLIB_UFC;
898    }
899
900    // Round up:
901    if ((rm == FPLIB_RN && (error == 3 ||
902                            (error == 2 && (int_mant & 1)))) ||
903        (((rm == FPLIB_RP && !sgn) || (rm == FPLIB_RM && sgn)) && error)) {
904        ++int_mant;
905        if (int_mant == 1ULL << FP64_MANT_BITS) {
906            // Rounded up from denormalized to normalized
907            biased_exp = 1;
908        }
909        if (int_mant == 2ULL << FP64_MANT_BITS) {
910            // Rounded up to next exponent
911            ++biased_exp;
912            int_mant >>= 1;
913        }
914    }
915
916    // Handle rounding to odd aka Von Neumann rounding:
917    if (error && rm == FPRounding_ODD)
918        int_mant |= 1;
919
920    // Handle overflow:
921    if (biased_exp >= (int)FP64_EXP_INF) {
922        *flags |= FPLIB_OFC | FPLIB_IXC;
923        if (rm == FPLIB_RN || (rm == FPLIB_RP && !sgn) ||
924            (rm == FPLIB_RM && sgn)) {
925            return fp64_infinity(sgn);
926        } else {
927            return fp64_max_normal(sgn);
928        }
929    }
930
931    if (error) {
932        *flags |= FPLIB_IXC;
933    }
934
935    return fp64_pack(sgn, biased_exp, int_mant);
936}
937
938static uint64_t
939fp64_round(int sgn, int exp, uint64_t mnt, int mode, int *flags)
940{
941    return fp64_round_(sgn, exp, mnt, mode & 3, mode, flags);
942}
943
944static int
945fp16_compare_eq(uint16_t a, uint16_t b, int mode, int *flags)
946{
947    int a_sgn, a_exp, b_sgn, b_exp;
948    uint16_t a_mnt, b_mnt;
949
950    fp16_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
951    fp16_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
952
953    if (fp16_is_NaN(a_exp, a_mnt) ||
954        fp16_is_NaN(b_exp, b_mnt)) {
955        if (fp16_is_signalling_NaN(a_exp, a_mnt) ||
956            fp16_is_signalling_NaN(b_exp, b_mnt))
957            *flags |= FPLIB_IOC;
958        return 0;
959    }
960    return a == b || (!a_mnt && !b_mnt);
961}
962
963static int
964fp16_compare_ge(uint16_t a, uint16_t b, int mode, int *flags)
965{
966    int a_sgn, a_exp, b_sgn, b_exp;
967    uint16_t a_mnt, b_mnt;
968
969    fp16_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
970    fp16_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
971
972    if (fp16_is_NaN(a_exp, a_mnt) ||
973        fp16_is_NaN(b_exp, b_mnt)) {
974        *flags |= FPLIB_IOC;
975        return 0;
976    }
977    if (!a_mnt && !b_mnt)
978        return 1;
979    if (a_sgn != b_sgn)
980        return b_sgn;
981    if (a_exp != b_exp)
982        return a_sgn ^ (a_exp > b_exp);
983    if (a_mnt != b_mnt)
984        return a_sgn ^ (a_mnt > b_mnt);
985    return 1;
986}
987
988static int
989fp16_compare_gt(uint16_t a, uint16_t b, int mode, int *flags)
990{
991    int a_sgn, a_exp, b_sgn, b_exp;
992    uint16_t a_mnt, b_mnt;
993
994    fp16_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
995    fp16_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
996
997    if (fp16_is_NaN(a_exp, a_mnt) ||
998        fp16_is_NaN(b_exp, b_mnt)) {
999        *flags |= FPLIB_IOC;
1000        return 0;
1001    }
1002    if (!a_mnt && !b_mnt)
1003        return 0;
1004    if (a_sgn != b_sgn)
1005        return b_sgn;
1006    if (a_exp != b_exp)
1007        return a_sgn ^ (a_exp > b_exp);
1008    if (a_mnt != b_mnt)
1009        return a_sgn ^ (a_mnt > b_mnt);
1010    return 0;
1011}
1012
1013static int
1014fp16_compare_un(uint16_t a, uint16_t b, int mode, int *flags)
1015{
1016    int a_sgn, a_exp, b_sgn, b_exp;
1017    uint16_t a_mnt, b_mnt;
1018
1019    fp16_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1020    fp16_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1021
1022    if (fp16_is_NaN(a_exp, a_mnt) ||
1023        fp16_is_NaN(b_exp, b_mnt)) {
1024        if (fp16_is_signalling_NaN(a_exp, a_mnt) ||
1025            fp16_is_signalling_NaN(b_exp, b_mnt))
1026            *flags |= FPLIB_IOC;
1027        return 1;
1028    }
1029    return 0;
1030}
1031
1032static int
1033fp32_compare_eq(uint32_t a, uint32_t b, int mode, int *flags)
1034{
1035    int a_sgn, a_exp, b_sgn, b_exp;
1036    uint32_t a_mnt, b_mnt;
1037
1038    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1039    fp32_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1040
1041    if (fp32_is_NaN(a_exp, a_mnt) ||
1042        fp32_is_NaN(b_exp, b_mnt)) {
1043        if (fp32_is_signalling_NaN(a_exp, a_mnt) ||
1044            fp32_is_signalling_NaN(b_exp, b_mnt))
1045            *flags |= FPLIB_IOC;
1046        return 0;
1047    }
1048    return a == b || (!a_mnt && !b_mnt);
1049}
1050
1051static int
1052fp32_compare_ge(uint32_t a, uint32_t b, int mode, int *flags)
1053{
1054    int a_sgn, a_exp, b_sgn, b_exp;
1055    uint32_t a_mnt, b_mnt;
1056
1057    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1058    fp32_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1059
1060    if (fp32_is_NaN(a_exp, a_mnt) ||
1061        fp32_is_NaN(b_exp, b_mnt)) {
1062        *flags |= FPLIB_IOC;
1063        return 0;
1064    }
1065    if (!a_mnt && !b_mnt)
1066        return 1;
1067    if (a_sgn != b_sgn)
1068        return b_sgn;
1069    if (a_exp != b_exp)
1070        return a_sgn ^ (a_exp > b_exp);
1071    if (a_mnt != b_mnt)
1072        return a_sgn ^ (a_mnt > b_mnt);
1073    return 1;
1074}
1075
1076static int
1077fp32_compare_gt(uint32_t a, uint32_t b, int mode, int *flags)
1078{
1079    int a_sgn, a_exp, b_sgn, b_exp;
1080    uint32_t a_mnt, b_mnt;
1081
1082    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1083    fp32_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1084
1085    if (fp32_is_NaN(a_exp, a_mnt) ||
1086        fp32_is_NaN(b_exp, b_mnt)) {
1087        *flags |= FPLIB_IOC;
1088        return 0;
1089    }
1090    if (!a_mnt && !b_mnt)
1091        return 0;
1092    if (a_sgn != b_sgn)
1093        return b_sgn;
1094    if (a_exp != b_exp)
1095        return a_sgn ^ (a_exp > b_exp);
1096    if (a_mnt != b_mnt)
1097        return a_sgn ^ (a_mnt > b_mnt);
1098    return 0;
1099}
1100
1101static int
1102fp32_compare_un(uint32_t a, uint32_t b, int mode, int *flags)
1103{
1104    int a_sgn, a_exp, b_sgn, b_exp;
1105    uint32_t a_mnt, b_mnt;
1106
1107    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1108    fp32_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1109
1110    if (fp32_is_NaN(a_exp, a_mnt) ||
1111        fp32_is_NaN(b_exp, b_mnt)) {
1112        if (fp32_is_signalling_NaN(a_exp, a_mnt) ||
1113            fp32_is_signalling_NaN(b_exp, b_mnt))
1114            *flags |= FPLIB_IOC;
1115        return 1;
1116    }
1117    return 0;
1118}
1119
1120static int
1121fp64_compare_eq(uint64_t a, uint64_t b, int mode, int *flags)
1122{
1123    int a_sgn, a_exp, b_sgn, b_exp;
1124    uint64_t a_mnt, b_mnt;
1125
1126    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1127    fp64_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1128
1129    if (fp64_is_NaN(a_exp, a_mnt) ||
1130        fp64_is_NaN(b_exp, b_mnt)) {
1131        if (fp64_is_signalling_NaN(a_exp, a_mnt) ||
1132            fp64_is_signalling_NaN(b_exp, b_mnt))
1133            *flags |= FPLIB_IOC;
1134        return 0;
1135    }
1136    return a == b || (!a_mnt && !b_mnt);
1137}
1138
1139static int
1140fp64_compare_ge(uint64_t a, uint64_t b, int mode, int *flags)
1141{
1142    int a_sgn, a_exp, b_sgn, b_exp;
1143    uint64_t a_mnt, b_mnt;
1144
1145    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1146    fp64_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1147
1148    if (fp64_is_NaN(a_exp, a_mnt) ||
1149        fp64_is_NaN(b_exp, b_mnt)) {
1150        *flags |= FPLIB_IOC;
1151        return 0;
1152    }
1153    if (!a_mnt && !b_mnt)
1154        return 1;
1155    if (a_sgn != b_sgn)
1156        return b_sgn;
1157    if (a_exp != b_exp)
1158        return a_sgn ^ (a_exp > b_exp);
1159    if (a_mnt != b_mnt)
1160        return a_sgn ^ (a_mnt > b_mnt);
1161    return 1;
1162}
1163
1164static int
1165fp64_compare_gt(uint64_t a, uint64_t b, int mode, int *flags)
1166{
1167    int a_sgn, a_exp, b_sgn, b_exp;
1168    uint64_t a_mnt, b_mnt;
1169
1170    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1171    fp64_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1172
1173    if (fp64_is_NaN(a_exp, a_mnt) ||
1174        fp64_is_NaN(b_exp, b_mnt)) {
1175        *flags |= FPLIB_IOC;
1176        return 0;
1177    }
1178    if (!a_mnt && !b_mnt)
1179        return 0;
1180    if (a_sgn != b_sgn)
1181        return b_sgn;
1182    if (a_exp != b_exp)
1183        return a_sgn ^ (a_exp > b_exp);
1184    if (a_mnt != b_mnt)
1185        return a_sgn ^ (a_mnt > b_mnt);
1186    return 0;
1187}
1188
1189static int
1190fp64_compare_un(uint64_t a, uint64_t b, int mode, int *flags)
1191{
1192    int a_sgn, a_exp, b_sgn, b_exp;
1193    uint64_t a_mnt, b_mnt;
1194
1195    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1196    fp64_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1197
1198    if (fp64_is_NaN(a_exp, a_mnt) ||
1199        fp64_is_NaN(b_exp, b_mnt)) {
1200        if (fp64_is_signalling_NaN(a_exp, a_mnt) ||
1201            fp64_is_signalling_NaN(b_exp, b_mnt))
1202            *flags |= FPLIB_IOC;
1203        return 1;
1204    }
1205    return 0;
1206}
1207
1208static uint16_t
1209fp16_add(uint16_t a, uint16_t b, int neg, int mode, int *flags)
1210{
1211    int a_sgn, a_exp, b_sgn, b_exp, x_sgn, x_exp;
1212    uint16_t a_mnt, b_mnt, x, x_mnt;
1213
1214    fp16_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1215    fp16_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1216
1217    if ((x = fp16_process_NaNs(a, b, mode, flags))) {
1218        return x;
1219    }
1220
1221    b_sgn ^= neg;
1222
1223    // Handle infinities and zeroes:
1224    if (a_exp == FP16_EXP_INF && b_exp == FP16_EXP_INF && a_sgn != b_sgn) {
1225        *flags |= FPLIB_IOC;
1226        return fp16_defaultNaN();
1227    } else if (a_exp == FP16_EXP_INF) {
1228        return fp16_infinity(a_sgn);
1229    } else if (b_exp == FP16_EXP_INF) {
1230        return fp16_infinity(b_sgn);
1231    } else if (!a_mnt && !b_mnt && a_sgn == b_sgn) {
1232        return fp16_zero(a_sgn);
1233    }
1234
1235    a_mnt <<= 3;
1236    b_mnt <<= 3;
1237    if (a_exp >= b_exp) {
1238        b_mnt = (lsr16(b_mnt, a_exp - b_exp) |
1239                 !!(b_mnt & (lsl16(1, a_exp - b_exp) - 1)));
1240        b_exp = a_exp;
1241    } else {
1242        a_mnt = (lsr16(a_mnt, b_exp - a_exp) |
1243                 !!(a_mnt & (lsl16(1, b_exp - a_exp) - 1)));
1244        a_exp = b_exp;
1245    }
1246    x_sgn = a_sgn;
1247    x_exp = a_exp;
1248    if (a_sgn == b_sgn) {
1249        x_mnt = a_mnt + b_mnt;
1250    } else if (a_mnt >= b_mnt) {
1251        x_mnt = a_mnt - b_mnt;
1252    } else {
1253        x_sgn ^= 1;
1254        x_mnt = b_mnt - a_mnt;
1255    }
1256
1257    if (!x_mnt) {
1258        // Sign of exact zero result depends on rounding mode
1259        return fp16_zero((mode & 3) == 2);
1260    }
1261
1262    x_mnt = fp16_normalise(x_mnt, &x_exp);
1263
1264    return fp16_round(x_sgn, x_exp + FP16_EXP_BITS - 3, x_mnt << 1,
1265                      mode, flags);
1266}
1267
1268static uint32_t
1269fp32_add(uint32_t a, uint32_t b, int neg, int mode, int *flags)
1270{
1271    int a_sgn, a_exp, b_sgn, b_exp, x_sgn, x_exp;
1272    uint32_t a_mnt, b_mnt, x, x_mnt;
1273
1274    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1275    fp32_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1276
1277    if ((x = fp32_process_NaNs(a, b, mode, flags))) {
1278        return x;
1279    }
1280
1281    b_sgn ^= neg;
1282
1283    // Handle infinities and zeroes:
1284    if (a_exp == FP32_EXP_INF && b_exp == FP32_EXP_INF && a_sgn != b_sgn) {
1285        *flags |= FPLIB_IOC;
1286        return fp32_defaultNaN();
1287    } else if (a_exp == FP32_EXP_INF) {
1288        return fp32_infinity(a_sgn);
1289    } else if (b_exp == FP32_EXP_INF) {
1290        return fp32_infinity(b_sgn);
1291    } else if (!a_mnt && !b_mnt && a_sgn == b_sgn) {
1292        return fp32_zero(a_sgn);
1293    }
1294
1295    a_mnt <<= 3;
1296    b_mnt <<= 3;
1297    if (a_exp >= b_exp) {
1298        b_mnt = (lsr32(b_mnt, a_exp - b_exp) |
1299                 !!(b_mnt & (lsl32(1, a_exp - b_exp) - 1)));
1300        b_exp = a_exp;
1301    } else {
1302        a_mnt = (lsr32(a_mnt, b_exp - a_exp) |
1303                 !!(a_mnt & (lsl32(1, b_exp - a_exp) - 1)));
1304        a_exp = b_exp;
1305    }
1306    x_sgn = a_sgn;
1307    x_exp = a_exp;
1308    if (a_sgn == b_sgn) {
1309        x_mnt = a_mnt + b_mnt;
1310    } else if (a_mnt >= b_mnt) {
1311        x_mnt = a_mnt - b_mnt;
1312    } else {
1313        x_sgn ^= 1;
1314        x_mnt = b_mnt - a_mnt;
1315    }
1316
1317    if (!x_mnt) {
1318        // Sign of exact zero result depends on rounding mode
1319        return fp32_zero((mode & 3) == 2);
1320    }
1321
1322    x_mnt = fp32_normalise(x_mnt, &x_exp);
1323
1324    return fp32_round(x_sgn, x_exp + FP32_EXP_BITS - 3, x_mnt << 1,
1325                      mode, flags);
1326}
1327
1328static uint64_t
1329fp64_add(uint64_t a, uint64_t b, int neg, int mode, int *flags)
1330{
1331    int a_sgn, a_exp, b_sgn, b_exp, x_sgn, x_exp;
1332    uint64_t a_mnt, b_mnt, x, x_mnt;
1333
1334    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1335    fp64_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1336
1337    if ((x = fp64_process_NaNs(a, b, mode, flags))) {
1338        return x;
1339    }
1340
1341    b_sgn ^= neg;
1342
1343    // Handle infinities and zeroes:
1344    if (a_exp == FP64_EXP_INF && b_exp == FP64_EXP_INF && a_sgn != b_sgn) {
1345        *flags |= FPLIB_IOC;
1346        return fp64_defaultNaN();
1347    } else if (a_exp == FP64_EXP_INF) {
1348        return fp64_infinity(a_sgn);
1349    } else if (b_exp == FP64_EXP_INF) {
1350        return fp64_infinity(b_sgn);
1351    } else if (!a_mnt && !b_mnt && a_sgn == b_sgn) {
1352        return fp64_zero(a_sgn);
1353    }
1354
1355    a_mnt <<= 3;
1356    b_mnt <<= 3;
1357    if (a_exp >= b_exp) {
1358        b_mnt = (lsr64(b_mnt, a_exp - b_exp) |
1359                 !!(b_mnt & (lsl64(1, a_exp - b_exp) - 1)));
1360        b_exp = a_exp;
1361    } else {
1362        a_mnt = (lsr64(a_mnt, b_exp - a_exp) |
1363                 !!(a_mnt & (lsl64(1, b_exp - a_exp) - 1)));
1364        a_exp = b_exp;
1365    }
1366    x_sgn = a_sgn;
1367    x_exp = a_exp;
1368    if (a_sgn == b_sgn) {
1369        x_mnt = a_mnt + b_mnt;
1370    } else if (a_mnt >= b_mnt) {
1371        x_mnt = a_mnt - b_mnt;
1372    } else {
1373        x_sgn ^= 1;
1374        x_mnt = b_mnt - a_mnt;
1375    }
1376
1377    if (!x_mnt) {
1378        // Sign of exact zero result depends on rounding mode
1379        return fp64_zero((mode & 3) == 2);
1380    }
1381
1382    x_mnt = fp64_normalise(x_mnt, &x_exp);
1383
1384    return fp64_round(x_sgn, x_exp + FP64_EXP_BITS - 3, x_mnt << 1,
1385                      mode, flags);
1386}
1387
1388static uint16_t
1389fp16_mul(uint16_t a, uint16_t b, int mode, int *flags)
1390{
1391    int a_sgn, a_exp, b_sgn, b_exp, x_sgn, x_exp;
1392    uint16_t a_mnt, b_mnt, x;
1393    uint32_t x_mnt;
1394
1395    fp16_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1396    fp16_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1397
1398    if ((x = fp16_process_NaNs(a, b, mode, flags))) {
1399        return x;
1400    }
1401
1402    // Handle infinities and zeroes:
1403    if ((a_exp == FP16_EXP_INF && !b_mnt) ||
1404        (b_exp == FP16_EXP_INF && !a_mnt)) {
1405        *flags |= FPLIB_IOC;
1406        return fp16_defaultNaN();
1407    } else if (a_exp == FP16_EXP_INF || b_exp == FP16_EXP_INF) {
1408        return fp16_infinity(a_sgn ^ b_sgn);
1409    } else if (!a_mnt || !b_mnt) {
1410        return fp16_zero(a_sgn ^ b_sgn);
1411    }
1412
1413    // Multiply and normalise:
1414    x_sgn = a_sgn ^ b_sgn;
1415    x_exp = a_exp + b_exp - FP16_EXP_BIAS + 2 * FP16_EXP_BITS + 1;
1416    x_mnt = (uint32_t)a_mnt * b_mnt;
1417    x_mnt = fp32_normalise(x_mnt, &x_exp);
1418
1419    // Convert to FP16_BITS bits, collapsing error into bottom bit:
1420    x_mnt = lsr32(x_mnt, FP16_BITS - 1) | !!lsl32(x_mnt, FP16_BITS + 1);
1421
1422    return fp16_round(x_sgn, x_exp, x_mnt, mode, flags);
1423}
1424
1425static uint32_t
1426fp32_mul(uint32_t a, uint32_t b, int mode, int *flags)
1427{
1428    int a_sgn, a_exp, b_sgn, b_exp, x_sgn, x_exp;
1429    uint32_t a_mnt, b_mnt, x;
1430    uint64_t x_mnt;
1431
1432    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1433    fp32_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1434
1435    if ((x = fp32_process_NaNs(a, b, mode, flags))) {
1436        return x;
1437    }
1438
1439    // Handle infinities and zeroes:
1440    if ((a_exp == FP32_EXP_INF && !b_mnt) ||
1441        (b_exp == FP32_EXP_INF && !a_mnt)) {
1442        *flags |= FPLIB_IOC;
1443        return fp32_defaultNaN();
1444    } else if (a_exp == FP32_EXP_INF || b_exp == FP32_EXP_INF) {
1445        return fp32_infinity(a_sgn ^ b_sgn);
1446    } else if (!a_mnt || !b_mnt) {
1447        return fp32_zero(a_sgn ^ b_sgn);
1448    }
1449
1450    // Multiply and normalise:
1451    x_sgn = a_sgn ^ b_sgn;
1452    x_exp = a_exp + b_exp - FP32_EXP_BIAS + 2 * FP32_EXP_BITS + 1;
1453    x_mnt = (uint64_t)a_mnt * b_mnt;
1454    x_mnt = fp64_normalise(x_mnt, &x_exp);
1455
1456    // Convert to FP32_BITS bits, collapsing error into bottom bit:
1457    x_mnt = lsr64(x_mnt, FP32_BITS - 1) | !!lsl64(x_mnt, FP32_BITS + 1);
1458
1459    return fp32_round(x_sgn, x_exp, x_mnt, mode, flags);
1460}
1461
1462static uint64_t
1463fp64_mul(uint64_t a, uint64_t b, int mode, int *flags)
1464{
1465    int a_sgn, a_exp, b_sgn, b_exp, x_sgn, x_exp;
1466    uint64_t a_mnt, b_mnt, x;
1467    uint64_t x0_mnt, x1_mnt;
1468
1469    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1470    fp64_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1471
1472    if ((x = fp64_process_NaNs(a, b, mode, flags))) {
1473        return x;
1474    }
1475
1476    // Handle infinities and zeroes:
1477    if ((a_exp == FP64_EXP_INF && !b_mnt) ||
1478        (b_exp == FP64_EXP_INF && !a_mnt)) {
1479        *flags |= FPLIB_IOC;
1480        return fp64_defaultNaN();
1481    } else if (a_exp == FP64_EXP_INF || b_exp == FP64_EXP_INF) {
1482        return fp64_infinity(a_sgn ^ b_sgn);
1483    } else if (!a_mnt || !b_mnt) {
1484        return fp64_zero(a_sgn ^ b_sgn);
1485    }
1486
1487    // Multiply and normalise:
1488    x_sgn = a_sgn ^ b_sgn;
1489    x_exp = a_exp + b_exp - FP64_EXP_BIAS + 2 * FP64_EXP_BITS + 1;
1490    mul62x62(&x0_mnt, &x1_mnt, a_mnt, b_mnt);
1491    fp128_normalise(&x0_mnt, &x1_mnt, &x_exp);
1492
1493    // Convert to FP64_BITS bits, collapsing error into bottom bit:
1494    x0_mnt = x1_mnt << 1 | !!x0_mnt;
1495
1496    return fp64_round(x_sgn, x_exp, x0_mnt, mode, flags);
1497}
1498
1499static uint16_t
1500fp16_muladd(uint16_t a, uint16_t b, uint16_t c, int scale,
1501            int mode, int *flags)
1502{
1503    int a_sgn, a_exp, b_sgn, b_exp, c_sgn, c_exp, x_sgn, x_exp, y_sgn, y_exp;
1504    uint16_t a_mnt, b_mnt, c_mnt, x;
1505    uint32_t x_mnt, y_mnt;
1506
1507    fp16_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1508    fp16_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1509    fp16_unpack(&c_sgn, &c_exp, &c_mnt, c, mode, flags);
1510
1511    x = fp16_process_NaNs3(a, b, c, mode, flags);
1512
1513    // Quiet NaN added to product of zero and infinity:
1514    if (fp16_is_quiet_NaN(a_exp, a_mnt) &&
1515        ((!b_mnt && fp16_is_infinity(c_exp, c_mnt)) ||
1516         (!c_mnt && fp16_is_infinity(b_exp, b_mnt)))) {
1517        x = fp16_defaultNaN();
1518        *flags |= FPLIB_IOC;
1519    }
1520
1521    if (x) {
1522        return x;
1523    }
1524
1525    // Handle infinities and zeroes:
1526    if ((b_exp == FP16_EXP_INF && !c_mnt) ||
1527        (c_exp == FP16_EXP_INF && !b_mnt) ||
1528        (a_exp == FP16_EXP_INF &&
1529         (b_exp == FP16_EXP_INF || c_exp == FP16_EXP_INF) &&
1530         (a_sgn != (b_sgn ^ c_sgn)))) {
1531        *flags |= FPLIB_IOC;
1532        return fp16_defaultNaN();
1533    }
1534    if (a_exp == FP16_EXP_INF)
1535        return fp16_infinity(a_sgn);
1536    if (b_exp == FP16_EXP_INF || c_exp == FP16_EXP_INF)
1537        return fp16_infinity(b_sgn ^ c_sgn);
1538    if (!a_mnt && (!b_mnt || !c_mnt) && a_sgn == (b_sgn ^ c_sgn))
1539        return fp16_zero(a_sgn);
1540
1541    x_sgn = a_sgn;
1542    x_exp = a_exp + 2 * FP16_EXP_BITS - 3;
1543    x_mnt = (uint32_t)a_mnt << (FP16_MANT_BITS + 4);
1544
1545    // Multiply:
1546    y_sgn = b_sgn ^ c_sgn;
1547    y_exp = b_exp + c_exp - FP16_EXP_BIAS + 2 * FP16_EXP_BITS + 1 - 3;
1548    y_mnt = (uint32_t)b_mnt * c_mnt << 3;
1549    if (!y_mnt) {
1550        y_exp = x_exp;
1551    }
1552
1553    // Add:
1554    if (x_exp >= y_exp) {
1555        y_mnt = (lsr32(y_mnt, x_exp - y_exp) |
1556                 !!(y_mnt & (lsl32(1, x_exp - y_exp) - 1)));
1557        y_exp = x_exp;
1558    } else {
1559        x_mnt = (lsr32(x_mnt, y_exp - x_exp) |
1560                 !!(x_mnt & (lsl32(1, y_exp - x_exp) - 1)));
1561        x_exp = y_exp;
1562    }
1563    if (x_sgn == y_sgn) {
1564        x_mnt = x_mnt + y_mnt;
1565    } else if (x_mnt >= y_mnt) {
1566        x_mnt = x_mnt - y_mnt;
1567    } else {
1568        x_sgn ^= 1;
1569        x_mnt = y_mnt - x_mnt;
1570    }
1571
1572    if (!x_mnt) {
1573        // Sign of exact zero result depends on rounding mode
1574        return fp16_zero((mode & 3) == 2);
1575    }
1576
1577    // Normalise into FP16_BITS bits, collapsing error into bottom bit:
1578    x_mnt = fp32_normalise(x_mnt, &x_exp);
1579    x_mnt = x_mnt >> (FP16_BITS - 1) | !!(uint16_t)(x_mnt << 1);
1580
1581    return fp16_round(x_sgn, x_exp + scale, x_mnt, mode, flags);
1582}
1583
1584static uint32_t
1585fp32_muladd(uint32_t a, uint32_t b, uint32_t c, int scale,
1586            int mode, int *flags)
1587{
1588    int a_sgn, a_exp, b_sgn, b_exp, c_sgn, c_exp, x_sgn, x_exp, y_sgn, y_exp;
1589    uint32_t a_mnt, b_mnt, c_mnt, x;
1590    uint64_t x_mnt, y_mnt;
1591
1592    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1593    fp32_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1594    fp32_unpack(&c_sgn, &c_exp, &c_mnt, c, mode, flags);
1595
1596    x = fp32_process_NaNs3(a, b, c, mode, flags);
1597
1598    // Quiet NaN added to product of zero and infinity:
1599    if (fp32_is_quiet_NaN(a_exp, a_mnt) &&
1600        ((!b_mnt && fp32_is_infinity(c_exp, c_mnt)) ||
1601         (!c_mnt && fp32_is_infinity(b_exp, b_mnt)))) {
1602        x = fp32_defaultNaN();
1603        *flags |= FPLIB_IOC;
1604    }
1605
1606    if (x) {
1607        return x;
1608    }
1609
1610    // Handle infinities and zeroes:
1611    if ((b_exp == FP32_EXP_INF && !c_mnt) ||
1612        (c_exp == FP32_EXP_INF && !b_mnt) ||
1613        (a_exp == FP32_EXP_INF &&
1614         (b_exp == FP32_EXP_INF || c_exp == FP32_EXP_INF) &&
1615         (a_sgn != (b_sgn ^ c_sgn)))) {
1616        *flags |= FPLIB_IOC;
1617        return fp32_defaultNaN();
1618    }
1619    if (a_exp == FP32_EXP_INF)
1620        return fp32_infinity(a_sgn);
1621    if (b_exp == FP32_EXP_INF || c_exp == FP32_EXP_INF)
1622        return fp32_infinity(b_sgn ^ c_sgn);
1623    if (!a_mnt && (!b_mnt || !c_mnt) && a_sgn == (b_sgn ^ c_sgn))
1624        return fp32_zero(a_sgn);
1625
1626    x_sgn = a_sgn;
1627    x_exp = a_exp + 2 * FP32_EXP_BITS - 3;
1628    x_mnt = (uint64_t)a_mnt << (FP32_MANT_BITS + 4);
1629
1630    // Multiply:
1631    y_sgn = b_sgn ^ c_sgn;
1632    y_exp = b_exp + c_exp - FP32_EXP_BIAS + 2 * FP32_EXP_BITS + 1 - 3;
1633    y_mnt = (uint64_t)b_mnt * c_mnt << 3;
1634    if (!y_mnt) {
1635        y_exp = x_exp;
1636    }
1637
1638    // Add:
1639    if (x_exp >= y_exp) {
1640        y_mnt = (lsr64(y_mnt, x_exp - y_exp) |
1641                 !!(y_mnt & (lsl64(1, x_exp - y_exp) - 1)));
1642        y_exp = x_exp;
1643    } else {
1644        x_mnt = (lsr64(x_mnt, y_exp - x_exp) |
1645                 !!(x_mnt & (lsl64(1, y_exp - x_exp) - 1)));
1646        x_exp = y_exp;
1647    }
1648    if (x_sgn == y_sgn) {
1649        x_mnt = x_mnt + y_mnt;
1650    } else if (x_mnt >= y_mnt) {
1651        x_mnt = x_mnt - y_mnt;
1652    } else {
1653        x_sgn ^= 1;
1654        x_mnt = y_mnt - x_mnt;
1655    }
1656
1657    if (!x_mnt) {
1658        // Sign of exact zero result depends on rounding mode
1659        return fp32_zero((mode & 3) == 2);
1660    }
1661
1662    // Normalise into FP32_BITS bits, collapsing error into bottom bit:
1663    x_mnt = fp64_normalise(x_mnt, &x_exp);
1664    x_mnt = x_mnt >> (FP32_BITS - 1) | !!(uint32_t)(x_mnt << 1);
1665
1666    return fp32_round(x_sgn, x_exp + scale, x_mnt, mode, flags);
1667}
1668
1669static uint64_t
1670fp64_muladd(uint64_t a, uint64_t b, uint64_t c, int scale,
1671            int mode, int *flags)
1672{
1673    int a_sgn, a_exp, b_sgn, b_exp, c_sgn, c_exp, x_sgn, x_exp, y_sgn, y_exp;
1674    uint64_t a_mnt, b_mnt, c_mnt, x;
1675    uint64_t x0_mnt, x1_mnt, y0_mnt, y1_mnt;
1676
1677    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1678    fp64_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1679    fp64_unpack(&c_sgn, &c_exp, &c_mnt, c, mode, flags);
1680
1681    x = fp64_process_NaNs3(a, b, c, mode, flags);
1682
1683    // Quiet NaN added to product of zero and infinity:
1684    if (fp64_is_quiet_NaN(a_exp, a_mnt) &&
1685        ((!b_mnt && fp64_is_infinity(c_exp, c_mnt)) ||
1686         (!c_mnt && fp64_is_infinity(b_exp, b_mnt)))) {
1687        x = fp64_defaultNaN();
1688        *flags |= FPLIB_IOC;
1689    }
1690
1691    if (x) {
1692        return x;
1693    }
1694
1695    // Handle infinities and zeroes:
1696    if ((b_exp == FP64_EXP_INF && !c_mnt) ||
1697        (c_exp == FP64_EXP_INF && !b_mnt) ||
1698        (a_exp == FP64_EXP_INF &&
1699         (b_exp == FP64_EXP_INF || c_exp == FP64_EXP_INF) &&
1700         (a_sgn != (b_sgn ^ c_sgn)))) {
1701        *flags |= FPLIB_IOC;
1702        return fp64_defaultNaN();
1703    }
1704    if (a_exp == FP64_EXP_INF)
1705        return fp64_infinity(a_sgn);
1706    if (b_exp == FP64_EXP_INF || c_exp == FP64_EXP_INF)
1707        return fp64_infinity(b_sgn ^ c_sgn);
1708    if (!a_mnt && (!b_mnt || !c_mnt) && a_sgn == (b_sgn ^ c_sgn))
1709        return fp64_zero(a_sgn);
1710
1711    x_sgn = a_sgn;
1712    x_exp = a_exp + FP64_EXP_BITS;
1713    x0_mnt = 0;
1714    x1_mnt = a_mnt;
1715
1716    // Multiply:
1717    y_sgn = b_sgn ^ c_sgn;
1718    y_exp = b_exp + c_exp - FP64_EXP_BIAS + 2 * FP64_EXP_BITS + 1 - 3;
1719    mul62x62(&y0_mnt, &y1_mnt, b_mnt, c_mnt << 3);
1720    if (!y0_mnt && !y1_mnt) {
1721        y_exp = x_exp;
1722    }
1723
1724    // Add:
1725    if (x_exp >= y_exp) {
1726        uint64_t t0, t1;
1727        lsl128(&t0, &t1, y0_mnt, y1_mnt,
1728               x_exp - y_exp < 128 ? 128 - (x_exp - y_exp) : 0);
1729        lsr128(&y0_mnt, &y1_mnt, y0_mnt, y1_mnt, x_exp - y_exp);
1730        y0_mnt |= !!(t0 | t1);
1731        y_exp = x_exp;
1732    } else {
1733        uint64_t t0, t1;
1734        lsl128(&t0, &t1, x0_mnt, x1_mnt,
1735               y_exp - x_exp < 128 ? 128 - (y_exp - x_exp) : 0);
1736        lsr128(&x0_mnt, &x1_mnt, x0_mnt, x1_mnt, y_exp - x_exp);
1737        x0_mnt |= !!(t0 | t1);
1738        x_exp = y_exp;
1739    }
1740    if (x_sgn == y_sgn) {
1741        add128(&x0_mnt, &x1_mnt, x0_mnt, x1_mnt, y0_mnt, y1_mnt);
1742    } else if (cmp128(x0_mnt, x1_mnt, y0_mnt, y1_mnt) >= 0) {
1743        sub128(&x0_mnt, &x1_mnt, x0_mnt, x1_mnt, y0_mnt, y1_mnt);
1744    } else {
1745        x_sgn ^= 1;
1746        sub128(&x0_mnt, &x1_mnt, y0_mnt, y1_mnt, x0_mnt, x1_mnt);
1747    }
1748
1749    if (!x0_mnt && !x1_mnt) {
1750        // Sign of exact zero result depends on rounding mode
1751        return fp64_zero((mode & 3) == 2);
1752    }
1753
1754    // Normalise into FP64_BITS bits, collapsing error into bottom bit:
1755    fp128_normalise(&x0_mnt, &x1_mnt, &x_exp);
1756    x0_mnt = x1_mnt << 1 | !!x0_mnt;
1757
1758    return fp64_round(x_sgn, x_exp + scale, x0_mnt, mode, flags);
1759}
1760
1761static uint16_t
1762fp16_div(uint16_t a, uint16_t b, int mode, int *flags)
1763{
1764    int a_sgn, a_exp, b_sgn, b_exp, x_sgn, x_exp;
1765    uint16_t a_mnt, b_mnt, x;
1766    uint32_t x_mnt;
1767
1768    fp16_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1769    fp16_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1770
1771    if ((x = fp16_process_NaNs(a, b, mode, flags)))
1772        return x;
1773
1774    // Handle infinities and zeroes:
1775    if ((a_exp == FP16_EXP_INF && b_exp == FP16_EXP_INF) ||
1776        (!a_mnt && !b_mnt)) {
1777        *flags |= FPLIB_IOC;
1778        return fp16_defaultNaN();
1779    }
1780    if (a_exp == FP16_EXP_INF || !b_mnt) {
1781        if (a_exp != FP16_EXP_INF)
1782            *flags |= FPLIB_DZC;
1783        return fp16_infinity(a_sgn ^ b_sgn);
1784    }
1785    if (!a_mnt || b_exp == FP16_EXP_INF)
1786        return fp16_zero(a_sgn ^ b_sgn);
1787
1788    // Divide, setting bottom bit if inexact:
1789    a_mnt = fp16_normalise(a_mnt, &a_exp);
1790    x_sgn = a_sgn ^ b_sgn;
1791    x_exp = a_exp - b_exp + (FP16_EXP_BIAS + FP16_BITS + 2 * FP16_EXP_BITS - 3);
1792    x_mnt = ((uint32_t)a_mnt << (FP16_MANT_BITS - FP16_EXP_BITS + 3)) / b_mnt;
1793    x_mnt |= (x_mnt * b_mnt !=
1794              (uint32_t)a_mnt << (FP16_MANT_BITS - FP16_EXP_BITS + 3));
1795
1796    // Normalise into FP16_BITS bits, collapsing error into bottom bit:
1797    x_mnt = fp32_normalise(x_mnt, &x_exp);
1798    x_mnt = x_mnt >> (FP16_BITS - 1) | !!(uint16_t)(x_mnt << 1);
1799
1800    return fp16_round(x_sgn, x_exp, x_mnt, mode, flags);
1801}
1802
1803static uint32_t
1804fp32_div(uint32_t a, uint32_t b, int mode, int *flags)
1805{
1806    int a_sgn, a_exp, b_sgn, b_exp, x_sgn, x_exp;
1807    uint32_t a_mnt, b_mnt, x;
1808    uint64_t x_mnt;
1809
1810    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1811    fp32_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1812
1813    if ((x = fp32_process_NaNs(a, b, mode, flags)))
1814        return x;
1815
1816    // Handle infinities and zeroes:
1817    if ((a_exp == FP32_EXP_INF && b_exp == FP32_EXP_INF) ||
1818        (!a_mnt && !b_mnt)) {
1819        *flags |= FPLIB_IOC;
1820        return fp32_defaultNaN();
1821    }
1822    if (a_exp == FP32_EXP_INF || !b_mnt) {
1823        if (a_exp != FP32_EXP_INF)
1824            *flags |= FPLIB_DZC;
1825        return fp32_infinity(a_sgn ^ b_sgn);
1826    }
1827    if (!a_mnt || b_exp == FP32_EXP_INF)
1828        return fp32_zero(a_sgn ^ b_sgn);
1829
1830    // Divide, setting bottom bit if inexact:
1831    a_mnt = fp32_normalise(a_mnt, &a_exp);
1832    x_sgn = a_sgn ^ b_sgn;
1833    x_exp = a_exp - b_exp + (FP32_EXP_BIAS + FP32_BITS + 2 * FP32_EXP_BITS - 3);
1834    x_mnt = ((uint64_t)a_mnt << (FP32_MANT_BITS - FP32_EXP_BITS + 3)) / b_mnt;
1835    x_mnt |= (x_mnt * b_mnt !=
1836              (uint64_t)a_mnt << (FP32_MANT_BITS - FP32_EXP_BITS + 3));
1837
1838    // Normalise into FP32_BITS bits, collapsing error into bottom bit:
1839    x_mnt = fp64_normalise(x_mnt, &x_exp);
1840    x_mnt = x_mnt >> (FP32_BITS - 1) | !!(uint32_t)(x_mnt << 1);
1841
1842    return fp32_round(x_sgn, x_exp, x_mnt, mode, flags);
1843}
1844
1845static uint64_t
1846fp64_div(uint64_t a, uint64_t b, int mode, int *flags)
1847{
1848    int a_sgn, a_exp, b_sgn, b_exp, x_sgn, x_exp, c;
1849    uint64_t a_mnt, b_mnt, x, x_mnt, x0_mnt, x1_mnt;
1850
1851    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1852    fp64_unpack(&b_sgn, &b_exp, &b_mnt, b, mode, flags);
1853
1854    if ((x = fp64_process_NaNs(a, b, mode, flags)))
1855        return x;
1856
1857    // Handle infinities and zeroes:
1858    if ((a_exp == FP64_EXP_INF && b_exp == FP64_EXP_INF) ||
1859        (!a_mnt && !b_mnt)) {
1860        *flags |= FPLIB_IOC;
1861        return fp64_defaultNaN();
1862    }
1863    if (a_exp == FP64_EXP_INF || !b_mnt) {
1864        if (a_exp != FP64_EXP_INF)
1865            *flags |= FPLIB_DZC;
1866        return fp64_infinity(a_sgn ^ b_sgn);
1867    }
1868    if (!a_mnt || b_exp == FP64_EXP_INF)
1869        return fp64_zero(a_sgn ^ b_sgn);
1870
1871    // Find reciprocal of divisor with Newton-Raphson:
1872    a_mnt = fp64_normalise(a_mnt, &a_exp);
1873    b_mnt = fp64_normalise(b_mnt, &b_exp);
1874    x_mnt = ~(uint64_t)0 / (b_mnt >> 31);
1875    mul64x32(&x0_mnt, &x1_mnt, b_mnt, x_mnt);
1876    sub128(&x0_mnt, &x1_mnt, 0, (uint64_t)1 << 32, x0_mnt, x1_mnt);
1877    lsr128(&x0_mnt, &x1_mnt, x0_mnt, x1_mnt, 32);
1878    mul64x32(&x0_mnt, &x1_mnt, x0_mnt, x_mnt);
1879    lsr128(&x0_mnt, &x1_mnt, x0_mnt, x1_mnt, 33);
1880
1881    // Multiply by dividend:
1882    x_sgn = a_sgn ^ b_sgn;
1883    x_exp = a_exp - b_exp + FP64_EXP_BIAS + 8;
1884    mul62x62(&x0_mnt, &x1_mnt, x0_mnt, a_mnt >> 2);
1885    lsr128(&x0_mnt, &x1_mnt, x0_mnt, x1_mnt, 4);
1886    x_mnt = x1_mnt;
1887
1888    // This is an underestimate, so try adding one:
1889    mul62x62(&x0_mnt, &x1_mnt, b_mnt >> 2, x_mnt + 1);
1890    c = cmp128(x0_mnt, x1_mnt, 0, a_mnt >> 11);
1891    if (c <= 0) {
1892        ++x_mnt;
1893    }
1894
1895    x_mnt = fp64_normalise(x_mnt, &x_exp);
1896
1897    return fp64_round(x_sgn, x_exp, x_mnt << 1 | !!c, mode, flags);
1898}
1899
1900static void
1901set_fpscr0(FPSCR &fpscr, int flags)
1902{
1903    if (flags & FPLIB_IDC) {
1904        fpscr.idc = 1;
1905    }
1906    if (flags & FPLIB_IOC) {
1907        fpscr.ioc = 1;
1908    }
1909    if (flags & FPLIB_DZC) {
1910        fpscr.dzc = 1;
1911    }
1912    if (flags & FPLIB_OFC) {
1913        fpscr.ofc = 1;
1914    }
1915    if (flags & FPLIB_UFC) {
1916        fpscr.ufc = 1;
1917    }
1918    if (flags & FPLIB_IXC) {
1919        fpscr.ixc = 1;
1920    }
1921}
1922
1923static uint16_t
1924fp16_scale(uint16_t a, int16_t b, int mode, int *flags)
1925{
1926    int a_sgn, a_exp;
1927    uint16_t a_mnt;
1928
1929    fp16_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1930
1931    // Handle NaNs:
1932    if (fp16_is_NaN(a_exp, a_mnt)) {
1933        return fp16_process_NaN(a, mode, flags);
1934    }
1935
1936    // Handle zeroes:
1937    if (!a_mnt) {
1938        return fp16_zero(a_sgn);
1939    }
1940
1941    // Handle infinities:
1942    if (a_exp == FP16_EXP_INF) {
1943        return fp16_infinity(a_sgn);
1944    }
1945
1946    b = b < -300 ? -300 : b;
1947    b = b >  300 ?  300 : b;
1948    a_exp += b;
1949    a_mnt <<= 3;
1950
1951    a_mnt = fp16_normalise(a_mnt, &a_exp);
1952
1953    return fp16_round(a_sgn, a_exp + FP16_EXP_BITS - 3, a_mnt << 1,
1954                      mode, flags);
1955}
1956
1957static uint32_t
1958fp32_scale(uint32_t a, int32_t b, int mode, int *flags)
1959{
1960    int a_sgn, a_exp;
1961    uint32_t a_mnt;
1962
1963    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1964
1965    // Handle NaNs:
1966    if (fp32_is_NaN(a_exp, a_mnt)) {
1967        return fp32_process_NaN(a, mode, flags);
1968    }
1969
1970    // Handle zeroes:
1971    if (!a_mnt) {
1972        return fp32_zero(a_sgn);
1973    }
1974
1975    // Handle infinities:
1976    if (a_exp == FP32_EXP_INF) {
1977        return fp32_infinity(a_sgn);
1978    }
1979
1980    b = b < -300 ? -300 : b;
1981    b = b >  300 ?  300 : b;
1982    a_exp += b;
1983    a_mnt <<= 3;
1984
1985    a_mnt = fp32_normalise(a_mnt, &a_exp);
1986
1987    return fp32_round(a_sgn, a_exp + FP32_EXP_BITS - 3, a_mnt << 1,
1988                      mode, flags);
1989}
1990
1991static uint64_t
1992fp64_scale(uint64_t a, int64_t b, int mode, int *flags)
1993{
1994    int a_sgn, a_exp;
1995    uint64_t a_mnt;
1996
1997    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
1998
1999    // Handle NaNs:
2000    if (fp64_is_NaN(a_exp, a_mnt)) {
2001        return fp64_process_NaN(a, mode, flags);
2002    }
2003
2004    // Handle zeroes:
2005    if (!a_mnt) {
2006        return fp64_zero(a_sgn);
2007    }
2008
2009    // Handle infinities:
2010    if (a_exp == FP64_EXP_INF) {
2011        return fp64_infinity(a_sgn);
2012    }
2013
2014    b = b < -3000 ? -3000 : b;
2015    b = b >  3000 ?  3000 : b;
2016    a_exp += b;
2017    a_mnt <<= 3;
2018
2019    a_mnt = fp64_normalise(a_mnt, &a_exp);
2020
2021    return fp64_round(a_sgn, a_exp + FP64_EXP_BITS - 3, a_mnt << 1,
2022                      mode, flags);
2023}
2024
2025static uint16_t
2026fp16_sqrt(uint16_t a, int mode, int *flags)
2027{
2028    int a_sgn, a_exp, x_sgn, x_exp;
2029    uint16_t a_mnt, x_mnt;
2030    uint32_t x, t0, t1;
2031
2032    fp16_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
2033
2034    // Handle NaNs:
2035    if (fp16_is_NaN(a_exp, a_mnt))
2036        return fp16_process_NaN(a, mode, flags);
2037
2038    // Handle infinities and zeroes:
2039    if (!a_mnt)
2040        return fp16_zero(a_sgn);
2041    if (a_exp == FP16_EXP_INF && !a_sgn)
2042        return fp16_infinity(a_sgn);
2043    if (a_sgn) {
2044        *flags |= FPLIB_IOC;
2045        return fp16_defaultNaN();
2046    }
2047
2048    a_mnt = fp16_normalise(a_mnt, &a_exp);
2049    if (a_exp & 1) {
2050        ++a_exp;
2051        a_mnt >>= 1;
2052    }
2053
2054    // x = (a * 3 + 5) / 8
2055    x = ((uint32_t)a_mnt << 14) + ((uint32_t)a_mnt << 13) + ((uint32_t)5 << 28);
2056
2057    // x = (a / x + x) / 2; // 8-bit accuracy
2058    x = (((uint32_t)a_mnt << 16) / (x >> 15) + (x >> 16)) << 15;
2059
2060    // x = (a / x + x) / 2; // 16-bit accuracy
2061    x = (((uint32_t)a_mnt << 16) / (x >> 15) + (x >> 16)) << 15;
2062
2063    x_sgn = 0;
2064    x_exp = (a_exp + 27) >> 1;
2065    x_mnt = ((x - (1 << 18)) >> 19) + 1;
2066    t1 = (uint32_t)x_mnt * x_mnt;
2067    t0 = (uint32_t)a_mnt << 9;
2068    if (t1 > t0) {
2069        --x_mnt;
2070    }
2071
2072    x_mnt = fp16_normalise(x_mnt, &x_exp);
2073
2074    return fp16_round(x_sgn, x_exp, x_mnt << 1 | (t1 != t0), mode, flags);
2075}
2076
2077static uint32_t
2078fp32_sqrt(uint32_t a, int mode, int *flags)
2079{
2080    int a_sgn, a_exp, x_sgn, x_exp;
2081    uint32_t a_mnt, x, x_mnt;
2082    uint64_t t0, t1;
2083
2084    fp32_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
2085
2086    // Handle NaNs:
2087    if (fp32_is_NaN(a_exp, a_mnt))
2088        return fp32_process_NaN(a, mode, flags);
2089
2090    // Handle infinities and zeroes:
2091    if (!a_mnt)
2092        return fp32_zero(a_sgn);
2093    if (a_exp == FP32_EXP_INF && !a_sgn)
2094        return fp32_infinity(a_sgn);
2095    if (a_sgn) {
2096        *flags |= FPLIB_IOC;
2097        return fp32_defaultNaN();
2098    }
2099
2100    a_mnt = fp32_normalise(a_mnt, &a_exp);
2101    if (!(a_exp & 1)) {
2102        ++a_exp;
2103        a_mnt >>= 1;
2104    }
2105
2106    // x = (a * 3 + 5) / 8
2107    x = (a_mnt >> 2) + (a_mnt >> 3) + ((uint32_t)5 << 28);
2108
2109    // x = (a / x + x) / 2; // 8-bit accuracy
2110    x = (a_mnt / (x >> 15) + (x >> 16)) << 15;
2111
2112    // x = (a / x + x) / 2; // 16-bit accuracy
2113    x = (a_mnt / (x >> 15) + (x >> 16)) << 15;
2114
2115    // x = (a / x + x) / 2; // 32-bit accuracy
2116    x = ((((uint64_t)a_mnt << 32) / x) >> 2) + (x >> 1);
2117
2118    x_sgn = 0;
2119    x_exp = (a_exp + 147) >> 1;
2120    x_mnt = ((x - (1 << 5)) >> 6) + 1;
2121    t1 = (uint64_t)x_mnt * x_mnt;
2122    t0 = (uint64_t)a_mnt << 19;
2123    if (t1 > t0) {
2124        --x_mnt;
2125    }
2126
2127    x_mnt = fp32_normalise(x_mnt, &x_exp);
2128
2129    return fp32_round(x_sgn, x_exp, x_mnt << 1 | (t1 != t0), mode, flags);
2130}
2131
2132static uint64_t
2133fp64_sqrt(uint64_t a, int mode, int *flags)
2134{
2135    int a_sgn, a_exp, x_sgn, x_exp, c;
2136    uint64_t a_mnt, x_mnt, r, x0, x1;
2137    uint32_t x;
2138
2139    fp64_unpack(&a_sgn, &a_exp, &a_mnt, a, mode, flags);
2140
2141    // Handle NaNs:
2142    if (fp64_is_NaN(a_exp, a_mnt))
2143        return fp64_process_NaN(a, mode, flags);
2144
2145    // Handle infinities and zeroes:
2146    if (!a_mnt)
2147        return fp64_zero(a_sgn);
2148    if (a_exp == FP64_EXP_INF && !a_sgn)
2149        return fp64_infinity(a_sgn);
2150    if (a_sgn) {
2151        *flags |= FPLIB_IOC;
2152        return fp64_defaultNaN();
2153    }
2154
2155    a_mnt = fp64_normalise(a_mnt, &a_exp);
2156    if (a_exp & 1) {
2157        ++a_exp;
2158        a_mnt >>= 1;
2159    }
2160
2161    // x = (a * 3 + 5) / 8
2162    x = (a_mnt >> 34) + (a_mnt >> 35) + ((uint32_t)5 << 28);
2163
2164    // x = (a / x + x) / 2; // 8-bit accuracy
2165    x = ((a_mnt >> 32) / (x >> 15) + (x >> 16)) << 15;
2166
2167    // x = (a / x + x) / 2; // 16-bit accuracy
2168    x = ((a_mnt >> 32) / (x >> 15) + (x >> 16)) << 15;
2169
2170    // x = (a / x + x) / 2; // 32-bit accuracy
2171    x = ((a_mnt / x) >> 2) + (x >> 1);
2172
2173    // r = 1 / x; // 32-bit accuracy
2174    r = ((uint64_t)1 << 62) / x;
2175
2176    // r = r * (2 - x * r); // 64-bit accuracy
2177    mul64x32(&x0, &x1, -(uint64_t)x * r << 1, r);
2178    lsr128(&x0, &x1, x0, x1, 31);
2179
2180    // x = (x + a * r) / 2; // 64-bit accuracy
2181    mul62x62(&x0, &x1, a_mnt >> 10, x0 >> 2);
2182    lsl128(&x0, &x1, x0, x1, 5);
2183    lsr128(&x0, &x1, x0, x1, 56);
2184
2185    x0 = ((uint64_t)x << 31) + (x0 >> 1);
2186
2187    x_sgn = 0;
2188    x_exp = (a_exp + 1053) >> 1;
2189    x_mnt = x0;
2190    x_mnt = ((x_mnt - (1 << 8)) >> 9) + 1;
2191    mul62x62(&x0, &x1, x_mnt, x_mnt);
2192    lsl128(&x0, &x1, x0, x1, 19);
2193    c = cmp128(x0, x1, 0, a_mnt);
2194    if (c > 0)
2195        --x_mnt;
2196
2197    x_mnt = fp64_normalise(x_mnt, &x_exp);
2198
2199    return fp64_round(x_sgn, x_exp, x_mnt << 1 | !!c, mode, flags);
2200}
2201
2202static int
2203modeConv(FPSCR fpscr)
2204{
2205    uint32_t x = (uint32_t)fpscr;
2206    return (x >> 22 & 0xf) | (x >> 19 & 1 ? FPLIB_FZ16 : 0);
2207    // AHP bit is ignored. Only fplibConvert uses AHP.
2208}
2209
2210static void
2211set_fpscr(FPSCR &fpscr, int flags)
2212{
2213    // translate back to FPSCR
2214    bool underflow = false;
2215    if (flags & FPLIB_IDC) {
2216        fpscr.idc = 1;
2217    }
2218    if (flags & FPLIB_IOC) {
2219        fpscr.ioc = 1;
2220    }
2221    if (flags & FPLIB_DZC) {
2222        fpscr.dzc = 1;
2223    }
2224    if (flags & FPLIB_OFC) {
2225        fpscr.ofc = 1;
2226    }
2227    if (flags & FPLIB_UFC) {
2228        underflow = true; //xx Why is this required?
2229        fpscr.ufc = 1;
2230    }
2231    if ((flags & FPLIB_IXC) && !(underflow && fpscr.fz)) {
2232        fpscr.ixc = 1;
2233    }
2234}
2235
2236template <>
2237bool
2238fplibCompareEQ(uint16_t a, uint16_t b, FPSCR &fpscr)
2239{
2240    int flags = 0;
2241    int x = fp16_compare_eq(a, b, modeConv(fpscr), &flags);
2242    set_fpscr(fpscr, flags);
2243    return x;
2244}
2245
2246template <>
2247bool
2248fplibCompareGE(uint16_t a, uint16_t b, FPSCR &fpscr)
2249{
2250    int flags = 0;
2251    int x = fp16_compare_ge(a, b, modeConv(fpscr), &flags);
2252    set_fpscr(fpscr, flags);
2253    return x;
2254}
2255
2256template <>
2257bool
2258fplibCompareGT(uint16_t a, uint16_t b, FPSCR &fpscr)
2259{
2260    int flags = 0;
2261    int x = fp16_compare_gt(a, b, modeConv(fpscr), &flags);
2262    set_fpscr(fpscr, flags);
2263    return x;
2264}
2265
2266template <>
2267bool
2268fplibCompareUN(uint16_t a, uint16_t b, FPSCR &fpscr)
2269{
2270    int flags = 0;
2271    int x = fp16_compare_un(a, b, modeConv(fpscr), &flags);
2272    set_fpscr(fpscr, flags);
2273    return x;
2274}
2275
2276template <>
2277bool
2278fplibCompareEQ(uint32_t a, uint32_t b, FPSCR &fpscr)
2279{
2280    int flags = 0;
2281    int x = fp32_compare_eq(a, b, modeConv(fpscr), &flags);
2282    set_fpscr(fpscr, flags);
2283    return x;
2284}
2285
2286template <>
2287bool
2288fplibCompareGE(uint32_t a, uint32_t b, FPSCR &fpscr)
2289{
2290    int flags = 0;
2291    int x = fp32_compare_ge(a, b, modeConv(fpscr), &flags);
2292    set_fpscr(fpscr, flags);
2293    return x;
2294}
2295
2296template <>
2297bool
2298fplibCompareGT(uint32_t a, uint32_t b, FPSCR &fpscr)
2299{
2300    int flags = 0;
2301    int x = fp32_compare_gt(a, b, modeConv(fpscr), &flags);
2302    set_fpscr(fpscr, flags);
2303    return x;
2304}
2305
2306template <>
2307bool
2308fplibCompareUN(uint32_t a, uint32_t b, FPSCR &fpscr)
2309{
2310    int flags = 0;
2311    int x = fp32_compare_un(a, b, modeConv(fpscr), &flags);
2312    set_fpscr(fpscr, flags);
2313    return x;
2314}
2315
2316template <>
2317bool
2318fplibCompareEQ(uint64_t a, uint64_t b, FPSCR &fpscr)
2319{
2320    int flags = 0;
2321    int x = fp64_compare_eq(a, b, modeConv(fpscr), &flags);
2322    set_fpscr(fpscr, flags);
2323    return x;
2324}
2325
2326template <>
2327bool
2328fplibCompareGE(uint64_t a, uint64_t b, FPSCR &fpscr)
2329{
2330    int flags = 0;
2331    int x = fp64_compare_ge(a, b, modeConv(fpscr), &flags);
2332    set_fpscr(fpscr, flags);
2333    return x;
2334}
2335
2336template <>
2337bool
2338fplibCompareGT(uint64_t a, uint64_t b, FPSCR &fpscr)
2339{
2340    int flags = 0;
2341    int x = fp64_compare_gt(a, b, modeConv(fpscr), &flags);
2342    set_fpscr(fpscr, flags);
2343    return x;
2344}
2345
2346template <>
2347bool
2348fplibCompareUN(uint64_t a, uint64_t b, FPSCR &fpscr)
2349{
2350    int flags = 0;
2351    int x = fp64_compare_un(a, b, modeConv(fpscr), &flags);
2352    set_fpscr(fpscr, flags);
2353    return x;
2354}
2355
2356template <>
2357uint16_t
2358fplibAbs(uint16_t op)
2359{
2360    return op & ~(1ULL << (FP16_BITS - 1));
2361}
2362
2363template <>
2364uint32_t
2365fplibAbs(uint32_t op)
2366{
2367    return op & ~(1ULL << (FP32_BITS - 1));
2368}
2369
2370template <>
2371uint64_t
2372fplibAbs(uint64_t op)
2373{
2374    return op & ~(1ULL << (FP64_BITS - 1));
2375}
2376
2377template <>
2378uint16_t
2379fplibAdd(uint16_t op1, uint16_t op2, FPSCR &fpscr)
2380{
2381    int flags = 0;
2382    uint16_t result = fp16_add(op1, op2, 0, modeConv(fpscr), &flags);
2383    set_fpscr0(fpscr, flags);
2384    return result;
2385}
2386
2387template <>
2388uint32_t
2389fplibAdd(uint32_t op1, uint32_t op2, FPSCR &fpscr)
2390{
2391    int flags = 0;
2392    uint32_t result = fp32_add(op1, op2, 0, modeConv(fpscr), &flags);
2393    set_fpscr0(fpscr, flags);
2394    return result;
2395}
2396
2397template <>
2398uint64_t
2399fplibAdd(uint64_t op1, uint64_t op2, FPSCR &fpscr)
2400{
2401    int flags = 0;
2402    uint64_t result = fp64_add(op1, op2, 0, modeConv(fpscr), &flags);
2403    set_fpscr0(fpscr, flags);
2404    return result;
2405}
2406
2407template <>
2408int
2409fplibCompare(uint16_t op1, uint16_t op2, bool signal_nans, FPSCR &fpscr)
2410{
2411    int mode = modeConv(fpscr);
2412    int flags = 0;
2413    int sgn1, exp1, sgn2, exp2, result;
2414    uint16_t mnt1, mnt2;
2415
2416    fp16_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
2417    fp16_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
2418
2419    if (fp16_is_NaN(exp1, mnt1) || fp16_is_NaN(exp2, mnt2)) {
2420        result = 3;
2421        if (fp16_is_signalling_NaN(exp1, mnt1) ||
2422            fp16_is_signalling_NaN(exp2, mnt2) || signal_nans)
2423            flags |= FPLIB_IOC;
2424    } else {
2425        if (op1 == op2 || (!mnt1 && !mnt2)) {
2426            result = 6;
2427        } else if (sgn1 != sgn2) {
2428            result = sgn1 ? 8 : 2;
2429        } else if (exp1 != exp2) {
2430            result = sgn1 ^ (exp1 < exp2) ? 8 : 2;
2431        } else {
2432            result = sgn1 ^ (mnt1 < mnt2) ? 8 : 2;
2433        }
2434    }
2435
2436    set_fpscr0(fpscr, flags);
2437
2438    return result;
2439}
2440
2441template <>
2442int
2443fplibCompare(uint32_t op1, uint32_t op2, bool signal_nans, FPSCR &fpscr)
2444{
2445    int mode = modeConv(fpscr);
2446    int flags = 0;
2447    int sgn1, exp1, sgn2, exp2, result;
2448    uint32_t mnt1, mnt2;
2449
2450    fp32_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
2451    fp32_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
2452
2453    if (fp32_is_NaN(exp1, mnt1) || fp32_is_NaN(exp2, mnt2)) {
2454        result = 3;
2455        if (fp32_is_signalling_NaN(exp1, mnt1) ||
2456            fp32_is_signalling_NaN(exp2, mnt2) || signal_nans)
2457            flags |= FPLIB_IOC;
2458    } else {
2459        if (op1 == op2 || (!mnt1 && !mnt2)) {
2460            result = 6;
2461        } else if (sgn1 != sgn2) {
2462            result = sgn1 ? 8 : 2;
2463        } else if (exp1 != exp2) {
2464            result = sgn1 ^ (exp1 < exp2) ? 8 : 2;
2465        } else {
2466            result = sgn1 ^ (mnt1 < mnt2) ? 8 : 2;
2467        }
2468    }
2469
2470    set_fpscr0(fpscr, flags);
2471
2472    return result;
2473}
2474
2475template <>
2476int
2477fplibCompare(uint64_t op1, uint64_t op2, bool signal_nans, FPSCR &fpscr)
2478{
2479    int mode = modeConv(fpscr);
2480    int flags = 0;
2481    int sgn1, exp1, sgn2, exp2, result;
2482    uint64_t mnt1, mnt2;
2483
2484    fp64_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
2485    fp64_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
2486
2487    if (fp64_is_NaN(exp1, mnt1) || fp64_is_NaN(exp2, mnt2)) {
2488        result = 3;
2489        if (fp64_is_signalling_NaN(exp1, mnt1) ||
2490            fp64_is_signalling_NaN(exp2, mnt2) || signal_nans)
2491            flags |= FPLIB_IOC;
2492    } else {
2493        if (op1 == op2 || (!mnt1 && !mnt2)) {
2494            result = 6;
2495        } else if (sgn1 != sgn2) {
2496            result = sgn1 ? 8 : 2;
2497        } else if (exp1 != exp2) {
2498            result = sgn1 ^ (exp1 < exp2) ? 8 : 2;
2499        } else {
2500            result = sgn1 ^ (mnt1 < mnt2) ? 8 : 2;
2501        }
2502    }
2503
2504    set_fpscr0(fpscr, flags);
2505
2506    return result;
2507}
2508
2509static uint16_t
2510fp16_FPConvertNaN_32(uint32_t op)
2511{
2512    return fp16_pack(op >> (FP32_BITS - 1), FP16_EXP_INF,
2513                     1ULL << (FP16_MANT_BITS - 1) |
2514                     op >> (FP32_MANT_BITS - FP16_MANT_BITS));
2515}
2516
2517static uint16_t
2518fp16_FPConvertNaN_64(uint64_t op)
2519{
2520    return fp16_pack(op >> (FP64_BITS - 1), FP16_EXP_INF,
2521                     1ULL << (FP16_MANT_BITS - 1) |
2522                     op >> (FP64_MANT_BITS - FP16_MANT_BITS));
2523}
2524
2525static uint32_t
2526fp32_FPConvertNaN_16(uint16_t op)
2527{
2528    return fp32_pack(op >> (FP16_BITS - 1), FP32_EXP_INF,
2529                     1ULL << (FP32_MANT_BITS - 1) |
2530                     (uint32_t)op << (FP32_MANT_BITS - FP16_MANT_BITS));
2531}
2532
2533static uint32_t
2534fp32_FPConvertNaN_64(uint64_t op)
2535{
2536    return fp32_pack(op >> (FP64_BITS - 1), FP32_EXP_INF,
2537                     1ULL << (FP32_MANT_BITS - 1) |
2538                     op >> (FP64_MANT_BITS - FP32_MANT_BITS));
2539}
2540
2541static uint64_t
2542fp64_FPConvertNaN_16(uint16_t op)
2543{
2544    return fp64_pack(op >> (FP16_BITS - 1), FP64_EXP_INF,
2545                     1ULL << (FP64_MANT_BITS - 1) |
2546                     (uint64_t)op << (FP64_MANT_BITS - FP16_MANT_BITS));
2547}
2548
2549static uint64_t
2550fp64_FPConvertNaN_32(uint32_t op)
2551{
2552    return fp64_pack(op >> (FP32_BITS - 1), FP64_EXP_INF,
2553                     1ULL << (FP64_MANT_BITS - 1) |
2554                     (uint64_t)op << (FP64_MANT_BITS - FP32_MANT_BITS));
2555}
2556
2557static uint16_t
2558fp16_FPOnePointFive(int sgn)
2559{
2560    return fp16_pack(sgn, FP16_EXP_BIAS, 1ULL << (FP16_MANT_BITS - 1));
2561}
2562
2563static uint32_t
2564fp32_FPOnePointFive(int sgn)
2565{
2566    return fp32_pack(sgn, FP32_EXP_BIAS, 1ULL << (FP32_MANT_BITS - 1));
2567}
2568
2569static uint64_t
2570fp64_FPOnePointFive(int sgn)
2571{
2572    return fp64_pack(sgn, FP64_EXP_BIAS, 1ULL << (FP64_MANT_BITS - 1));
2573}
2574
2575static uint16_t
2576fp16_FPThree(int sgn)
2577{
2578    return fp16_pack(sgn, FP16_EXP_BIAS + 1, 1ULL << (FP16_MANT_BITS - 1));
2579}
2580
2581static uint32_t
2582fp32_FPThree(int sgn)
2583{
2584    return fp32_pack(sgn, FP32_EXP_BIAS + 1, 1ULL << (FP32_MANT_BITS - 1));
2585}
2586
2587static uint64_t
2588fp64_FPThree(int sgn)
2589{
2590    return fp64_pack(sgn, FP64_EXP_BIAS + 1, 1ULL << (FP64_MANT_BITS - 1));
2591}
2592
2593static uint16_t
2594fp16_FPTwo(int sgn)
2595{
2596    return fp16_pack(sgn, FP16_EXP_BIAS + 1, 0);
2597}
2598
2599static uint32_t
2600fp32_FPTwo(int sgn)
2601{
2602    return fp32_pack(sgn, FP32_EXP_BIAS + 1, 0);
2603}
2604
2605static uint64_t
2606fp64_FPTwo(int sgn)
2607{
2608    return fp64_pack(sgn, FP64_EXP_BIAS + 1, 0);
2609}
2610
2611template <>
2612uint16_t
2613fplibConvert(uint32_t op, FPRounding rounding, FPSCR &fpscr)
2614{
2615    int mode = modeConv(fpscr);
2616    int flags = 0;
2617    int sgn, exp;
2618    uint32_t mnt;
2619    uint16_t result;
2620
2621    // Unpack floating-point operand optionally with flush-to-zero:
2622    fp32_unpack(&sgn, &exp, &mnt, op, mode, &flags);
2623
2624    bool alt_hp = fpscr.ahp;
2625
2626    if (fp32_is_NaN(exp, mnt)) {
2627        if (alt_hp) {
2628            result = fp16_zero(sgn);
2629        } else if (fpscr.dn) {
2630            result = fp16_defaultNaN();
2631        } else {
2632            result = fp16_FPConvertNaN_32(op);
2633        }
2634        if (!(mnt >> (FP32_MANT_BITS - 1) & 1) || alt_hp) {
2635            flags |= FPLIB_IOC;
2636        }
2637    } else if (exp == FP32_EXP_INF) {
2638        if (alt_hp) {
2639            result = ((uint16_t)sgn << (FP16_BITS - 1) |
2640                      ((1ULL << (FP16_BITS - 1)) - 1));
2641            flags |= FPLIB_IOC;
2642        } else {
2643            result = fp16_infinity(sgn);
2644        }
2645    } else if (!mnt) {
2646        result = fp16_zero(sgn);
2647    } else {
2648        result =
2649            fp16_round_(sgn, exp - FP32_EXP_BIAS + FP16_EXP_BIAS,
2650                        mnt >> (FP32_MANT_BITS - FP16_BITS) |
2651                        !!(mnt & ((1ULL << (FP32_MANT_BITS - FP16_BITS)) - 1)),
2652                        rounding, (mode & 0xf) | alt_hp << 4, &flags);
2653    }
2654
2655    set_fpscr0(fpscr, flags);
2656
2657    return result;
2658}
2659
2660template <>
2661uint16_t
2662fplibConvert(uint64_t op, FPRounding rounding, FPSCR &fpscr)
2663{
2664    int mode = modeConv(fpscr);
2665    int flags = 0;
2666    int sgn, exp;
2667    uint64_t mnt;
2668    uint16_t result;
2669
2670    // Unpack floating-point operand optionally with flush-to-zero:
2671    fp64_unpack(&sgn, &exp, &mnt, op, mode, &flags);
2672
2673    bool alt_hp = fpscr.ahp;
2674
2675    if (fp64_is_NaN(exp, mnt)) {
2676        if (alt_hp) {
2677            result = fp16_zero(sgn);
2678        } else if (fpscr.dn) {
2679            result = fp16_defaultNaN();
2680        } else {
2681            result = fp16_FPConvertNaN_64(op);
2682        }
2683        if (!(mnt >> (FP64_MANT_BITS - 1) & 1) || alt_hp) {
2684            flags |= FPLIB_IOC;
2685        }
2686    } else if (exp == FP64_EXP_INF) {
2687        if (alt_hp) {
2688            result = ((uint16_t)sgn << (FP16_BITS - 1) |
2689                      ((1ULL << (FP16_BITS - 1)) - 1));
2690            flags |= FPLIB_IOC;
2691        } else {
2692            result = fp16_infinity(sgn);
2693        }
2694    } else if (!mnt) {
2695        result = fp16_zero(sgn);
2696    } else {
2697        result =
2698            fp16_round_(sgn, exp - FP64_EXP_BIAS + FP16_EXP_BIAS,
2699                        mnt >> (FP64_MANT_BITS - FP16_BITS) |
2700                        !!(mnt & ((1ULL << (FP64_MANT_BITS - FP16_BITS)) - 1)),
2701                        rounding, (mode & 0xf) | alt_hp << 4, &flags);
2702    }
2703
2704    set_fpscr0(fpscr, flags);
2705
2706    return result;
2707}
2708
2709template <>
2710uint32_t
2711fplibConvert(uint16_t op, FPRounding rounding, FPSCR &fpscr)
2712{
2713    int mode = modeConv(fpscr);
2714    int flags = 0;
2715    int sgn, exp;
2716    uint16_t mnt;
2717    uint32_t result;
2718
2719    // Unpack floating-point operand optionally with flush-to-zero:
2720    fp16_unpack(&sgn, &exp, &mnt, op, mode & 0xf, &flags);
2721
2722    if (fp16_is_NaN(exp, mnt) && !fpscr.ahp) {
2723        if (fpscr.dn) {
2724            result = fp32_defaultNaN();
2725        } else {
2726            result = fp32_FPConvertNaN_16(op);
2727        }
2728        if (!(mnt >> (FP16_MANT_BITS - 1) & 1)) {
2729            flags |= FPLIB_IOC;
2730        }
2731    } else if (exp == FP16_EXP_INF && !fpscr.ahp) {
2732        result = fp32_infinity(sgn);
2733    } else if (!mnt) {
2734        result = fp32_zero(sgn);
2735    } else {
2736        mnt = fp16_normalise(mnt, &exp);
2737        result = fp32_pack(sgn, (exp - FP16_EXP_BIAS +
2738                                 FP32_EXP_BIAS + FP16_EXP_BITS),
2739                           (uint32_t)mnt << (FP32_MANT_BITS - FP16_BITS + 1));
2740    }
2741
2742    set_fpscr0(fpscr, flags);
2743
2744    return result;
2745}
2746
2747template <>
2748uint32_t
2749fplibConvert(uint64_t op, FPRounding rounding, FPSCR &fpscr)
2750{
2751    int mode = modeConv(fpscr);
2752    int flags = 0;
2753    int sgn, exp;
2754    uint64_t mnt;
2755    uint32_t result;
2756
2757    // Unpack floating-point operand optionally with flush-to-zero:
2758    fp64_unpack(&sgn, &exp, &mnt, op, mode, &flags);
2759
2760    if (fp64_is_NaN(exp, mnt)) {
2761        if (fpscr.dn) {
2762            result = fp32_defaultNaN();
2763        } else {
2764            result = fp32_FPConvertNaN_64(op);
2765        }
2766        if (!(mnt >> (FP64_MANT_BITS - 1) & 1)) {
2767            flags |= FPLIB_IOC;
2768        }
2769    } else if (exp == FP64_EXP_INF) {
2770        result = fp32_infinity(sgn);
2771    } else if (!mnt) {
2772        result = fp32_zero(sgn);
2773    } else {
2774        result =
2775            fp32_round_(sgn, exp - FP64_EXP_BIAS + FP32_EXP_BIAS,
2776                        mnt >> (FP64_MANT_BITS - FP32_BITS) |
2777                        !!(mnt & ((1ULL << (FP64_MANT_BITS - FP32_BITS)) - 1)),
2778                        rounding, mode, &flags);
2779    }
2780
2781    set_fpscr0(fpscr, flags);
2782
2783    return result;
2784}
2785
2786template <>
2787uint64_t
2788fplibConvert(uint16_t op, FPRounding rounding, FPSCR &fpscr)
2789{
2790    int mode = modeConv(fpscr);
2791    int flags = 0;
2792    int sgn, exp;
2793    uint16_t mnt;
2794    uint64_t result;
2795
2796    // Unpack floating-point operand optionally with flush-to-zero:
2797    fp16_unpack(&sgn, &exp, &mnt, op, mode & 0xf, &flags);
2798
2799    if (fp16_is_NaN(exp, mnt) && !fpscr.ahp) {
2800        if (fpscr.dn) {
2801            result = fp64_defaultNaN();
2802        } else {
2803            result = fp64_FPConvertNaN_16(op);
2804        }
2805        if (!(mnt >> (FP16_MANT_BITS - 1) & 1)) {
2806            flags |= FPLIB_IOC;
2807        }
2808    } else if (exp == FP16_EXP_INF && !fpscr.ahp) {
2809        result = fp64_infinity(sgn);
2810    } else if (!mnt) {
2811        result = fp64_zero(sgn);
2812    } else {
2813        mnt = fp16_normalise(mnt, &exp);
2814        result = fp64_pack(sgn, (exp - FP16_EXP_BIAS +
2815                                 FP64_EXP_BIAS + FP16_EXP_BITS),
2816                           (uint64_t)mnt << (FP64_MANT_BITS - FP16_BITS + 1));
2817    }
2818
2819    set_fpscr0(fpscr, flags);
2820
2821    return result;
2822}
2823
2824template <>
2825uint64_t
2826fplibConvert(uint32_t op, FPRounding rounding, FPSCR &fpscr)
2827{
2828    int mode = modeConv(fpscr);
2829    int flags = 0;
2830    int sgn, exp;
2831    uint32_t mnt;
2832    uint64_t result;
2833
2834    // Unpack floating-point operand optionally with flush-to-zero:
2835    fp32_unpack(&sgn, &exp, &mnt, op, mode, &flags);
2836
2837    if (fp32_is_NaN(exp, mnt)) {
2838        if (fpscr.dn) {
2839            result = fp64_defaultNaN();
2840        } else {
2841            result = fp64_FPConvertNaN_32(op);
2842        }
2843        if (!(mnt >> (FP32_MANT_BITS - 1) & 1)) {
2844            flags |= FPLIB_IOC;
2845        }
2846    } else if (exp == FP32_EXP_INF) {
2847        result = fp64_infinity(sgn);
2848    } else if (!mnt) {
2849        result = fp64_zero(sgn);
2850    } else {
2851        mnt = fp32_normalise(mnt, &exp);
2852        result = fp64_pack(sgn, (exp - FP32_EXP_BIAS +
2853                                 FP64_EXP_BIAS + FP32_EXP_BITS),
2854                           (uint64_t)mnt << (FP64_MANT_BITS - FP32_BITS + 1));
2855    }
2856
2857    set_fpscr0(fpscr, flags);
2858
2859    return result;
2860}
2861
2862template <>
2863uint16_t
2864fplibMulAdd(uint16_t addend, uint16_t op1, uint16_t op2, FPSCR &fpscr)
2865{
2866    int flags = 0;
2867    uint16_t result = fp16_muladd(addend, op1, op2, 0, modeConv(fpscr), &flags);
2868    set_fpscr0(fpscr, flags);
2869    return result;
2870}
2871
2872template <>
2873uint32_t
2874fplibMulAdd(uint32_t addend, uint32_t op1, uint32_t op2, FPSCR &fpscr)
2875{
2876    int flags = 0;
2877    uint32_t result = fp32_muladd(addend, op1, op2, 0, modeConv(fpscr), &flags);
2878    set_fpscr0(fpscr, flags);
2879    return result;
2880}
2881
2882template <>
2883uint64_t
2884fplibMulAdd(uint64_t addend, uint64_t op1, uint64_t op2, FPSCR &fpscr)
2885{
2886    int flags = 0;
2887    uint64_t result = fp64_muladd(addend, op1, op2, 0, modeConv(fpscr), &flags);
2888    set_fpscr0(fpscr, flags);
2889    return result;
2890}
2891
2892template <>
2893uint16_t
2894fplibDiv(uint16_t op1, uint16_t op2, FPSCR &fpscr)
2895{
2896    int flags = 0;
2897    uint16_t result = fp16_div(op1, op2, modeConv(fpscr), &flags);
2898    set_fpscr0(fpscr, flags);
2899    return result;
2900}
2901
2902template <>
2903uint32_t
2904fplibDiv(uint32_t op1, uint32_t op2, FPSCR &fpscr)
2905{
2906    int flags = 0;
2907    uint32_t result = fp32_div(op1, op2, modeConv(fpscr), &flags);
2908    set_fpscr0(fpscr, flags);
2909    return result;
2910}
2911
2912template <>
2913uint64_t
2914fplibDiv(uint64_t op1, uint64_t op2, FPSCR &fpscr)
2915{
2916    int flags = 0;
2917    uint64_t result = fp64_div(op1, op2, modeConv(fpscr), &flags);
2918    set_fpscr0(fpscr, flags);
2919    return result;
2920}
2921
2922template <>
2923uint16_t
2924fplibExpA(uint16_t op)
2925{
2926    static uint16_t coeff[32] = {
2927        0x0000,
2928        0x0016,
2929        0x002d,
2930        0x0045,
2931        0x005d,
2932        0x0075,
2933        0x008e,
2934        0x00a8,
2935        0x00c2,
2936        0x00dc,
2937        0x00f8,
2938        0x0114,
2939        0x0130,
2940        0x014d,
2941        0x016b,
2942        0x0189,
2943        0x01a8,
2944        0x01c8,
2945        0x01e8,
2946        0x0209,
2947        0x022b,
2948        0x024e,
2949        0x0271,
2950        0x0295,
2951        0x02ba,
2952        0x02e0,
2953        0x0306,
2954        0x032e,
2955        0x0356,
2956        0x037f,
2957        0x03a9,
2958        0x03d4
2959    };
2960    return ((((op >> 5) & ((1 << FP16_EXP_BITS) - 1)) << FP16_MANT_BITS) |
2961            coeff[op & ((1 << 5) - 1)]);
2962}
2963
2964template <>
2965uint32_t
2966fplibExpA(uint32_t op)
2967{
2968    static uint32_t coeff[64] = {
2969        0x000000,
2970        0x0164d2,
2971        0x02cd87,
2972        0x043a29,
2973        0x05aac3,
2974        0x071f62,
2975        0x08980f,
2976        0x0a14d5,
2977        0x0b95c2,
2978        0x0d1adf,
2979        0x0ea43a,
2980        0x1031dc,
2981        0x11c3d3,
2982        0x135a2b,
2983        0x14f4f0,
2984        0x16942d,
2985        0x1837f0,
2986        0x19e046,
2987        0x1b8d3a,
2988        0x1d3eda,
2989        0x1ef532,
2990        0x20b051,
2991        0x227043,
2992        0x243516,
2993        0x25fed7,
2994        0x27cd94,
2995        0x29a15b,
2996        0x2b7a3a,
2997        0x2d583f,
2998        0x2f3b79,
2999        0x3123f6,
3000        0x3311c4,
3001        0x3504f3,
3002        0x36fd92,
3003        0x38fbaf,
3004        0x3aff5b,
3005        0x3d08a4,
3006        0x3f179a,
3007        0x412c4d,
3008        0x4346cd,
3009        0x45672a,
3010        0x478d75,
3011        0x49b9be,
3012        0x4bec15,
3013        0x4e248c,
3014        0x506334,
3015        0x52a81e,
3016        0x54f35b,
3017        0x5744fd,
3018        0x599d16,
3019        0x5bfbb8,
3020        0x5e60f5,
3021        0x60ccdf,
3022        0x633f89,
3023        0x65b907,
3024        0x68396a,
3025        0x6ac0c7,
3026        0x6d4f30,
3027        0x6fe4ba,
3028        0x728177,
3029        0x75257d,
3030        0x77d0df,
3031        0x7a83b3,
3032        0x7d3e0c
3033    };
3034    return ((((op >> 6) & ((1 << FP32_EXP_BITS) - 1)) << FP32_MANT_BITS) |
3035            coeff[op & ((1 << 6) - 1)]);
3036}
3037
3038template <>
3039uint64_t
3040fplibExpA(uint64_t op)
3041{
3042    static uint64_t coeff[64] = {
3043        0x0000000000000ULL,
3044        0x02c9a3e778061ULL,
3045        0x059b0d3158574ULL,
3046        0x0874518759bc8ULL,
3047        0x0b5586cf9890fULL,
3048        0x0e3ec32d3d1a2ULL,
3049        0x11301d0125b51ULL,
3050        0x1429aaea92de0ULL,
3051        0x172b83c7d517bULL,
3052        0x1a35beb6fcb75ULL,
3053        0x1d4873168b9aaULL,
3054        0x2063b88628cd6ULL,
3055        0x2387a6e756238ULL,
3056        0x26b4565e27cddULL,
3057        0x29e9df51fdee1ULL,
3058        0x2d285a6e4030bULL,
3059        0x306fe0a31b715ULL,
3060        0x33c08b26416ffULL,
3061        0x371a7373aa9cbULL,
3062        0x3a7db34e59ff7ULL,
3063        0x3dea64c123422ULL,
3064        0x4160a21f72e2aULL,
3065        0x44e086061892dULL,
3066        0x486a2b5c13cd0ULL,
3067        0x4bfdad5362a27ULL,
3068        0x4f9b2769d2ca7ULL,
3069        0x5342b569d4f82ULL,
3070        0x56f4736b527daULL,
3071        0x5ab07dd485429ULL,
3072        0x5e76f15ad2148ULL,
3073        0x6247eb03a5585ULL,
3074        0x6623882552225ULL,
3075        0x6a09e667f3bcdULL,
3076        0x6dfb23c651a2fULL,
3077        0x71f75e8ec5f74ULL,
3078        0x75feb564267c9ULL,
3079        0x7a11473eb0187ULL,
3080        0x7e2f336cf4e62ULL,
3081        0x82589994cce13ULL,
3082        0x868d99b4492edULL,
3083        0x8ace5422aa0dbULL,
3084        0x8f1ae99157736ULL,
3085        0x93737b0cdc5e5ULL,
3086        0x97d829fde4e50ULL,
3087        0x9c49182a3f090ULL,
3088        0xa0c667b5de565ULL,
3089        0xa5503b23e255dULL,
3090        0xa9e6b5579fdbfULL,
3091        0xae89f995ad3adULL,
3092        0xb33a2b84f15fbULL,
3093        0xb7f76f2fb5e47ULL,
3094        0xbcc1e904bc1d2ULL,
3095        0xc199bdd85529cULL,
3096        0xc67f12e57d14bULL,
3097        0xcb720dcef9069ULL,
3098        0xd072d4a07897cULL,
3099        0xd5818dcfba487ULL,
3100        0xda9e603db3285ULL,
3101        0xdfc97337b9b5fULL,
3102        0xe502ee78b3ff6ULL,
3103        0xea4afa2a490daULL,
3104        0xefa1bee615a27ULL,
3105        0xf50765b6e4540ULL,
3106        0xfa7c1819e90d8ULL
3107    };
3108    return ((((op >> 6) & ((1 << FP64_EXP_BITS) - 1)) << FP64_MANT_BITS) |
3109            coeff[op & ((1 << 6) - 1)]);
3110}
3111
3112static uint16_t
3113fp16_repack(int sgn, int exp, uint16_t mnt)
3114{
3115    return fp16_pack(sgn, mnt >> FP16_MANT_BITS ? exp : 0, mnt);
3116}
3117
3118static uint32_t
3119fp32_repack(int sgn, int exp, uint32_t mnt)
3120{
3121    return fp32_pack(sgn, mnt >> FP32_MANT_BITS ? exp : 0, mnt);
3122}
3123
3124static uint64_t
3125fp64_repack(int sgn, int exp, uint64_t mnt)
3126{
3127    return fp64_pack(sgn, mnt >> FP64_MANT_BITS ? exp : 0, mnt);
3128}
3129
3130static void
3131fp16_minmaxnum(uint16_t *op1, uint16_t *op2, int sgn)
3132{
3133    // Treat a single quiet-NaN as +Infinity/-Infinity
3134    if (!((uint16_t)~(*op1 << 1) >> FP16_MANT_BITS) &&
3135        (uint16_t)~(*op2 << 1) >> FP16_MANT_BITS)
3136        *op1 = fp16_infinity(sgn);
3137    if (!((uint16_t)~(*op2 << 1) >> FP16_MANT_BITS) &&
3138        (uint16_t)~(*op1 << 1) >> FP16_MANT_BITS)
3139        *op2 = fp16_infinity(sgn);
3140}
3141
3142static void
3143fp32_minmaxnum(uint32_t *op1, uint32_t *op2, int sgn)
3144{
3145    // Treat a single quiet-NaN as +Infinity/-Infinity
3146    if (!((uint32_t)~(*op1 << 1) >> FP32_MANT_BITS) &&
3147        (uint32_t)~(*op2 << 1) >> FP32_MANT_BITS)
3148        *op1 = fp32_infinity(sgn);
3149    if (!((uint32_t)~(*op2 << 1) >> FP32_MANT_BITS) &&
3150        (uint32_t)~(*op1 << 1) >> FP32_MANT_BITS)
3151        *op2 = fp32_infinity(sgn);
3152}
3153
3154static void
3155fp64_minmaxnum(uint64_t *op1, uint64_t *op2, int sgn)
3156{
3157    // Treat a single quiet-NaN as +Infinity/-Infinity
3158    if (!((uint64_t)~(*op1 << 1) >> FP64_MANT_BITS) &&
3159        (uint64_t)~(*op2 << 1) >> FP64_MANT_BITS)
3160        *op1 = fp64_infinity(sgn);
3161    if (!((uint64_t)~(*op2 << 1) >> FP64_MANT_BITS) &&
3162        (uint64_t)~(*op1 << 1) >> FP64_MANT_BITS)
3163        *op2 = fp64_infinity(sgn);
3164}
3165
3166template <>
3167uint16_t
3168fplibMax(uint16_t op1, uint16_t op2, FPSCR &fpscr)
3169{
3170    int mode = modeConv(fpscr);
3171    int flags = 0;
3172    int sgn1, exp1, sgn2, exp2;
3173    uint16_t mnt1, mnt2, x, result;
3174
3175    fp16_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
3176    fp16_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
3177
3178    if ((x = fp16_process_NaNs(op1, op2, mode, &flags))) {
3179        result = x;
3180    } else {
3181        result = ((sgn1 != sgn2 ? sgn2 : sgn1 ^ (op1 > op2)) ?
3182                  fp16_repack(sgn1, exp1, mnt1) :
3183                  fp16_repack(sgn2, exp2, mnt2));
3184    }
3185    set_fpscr0(fpscr, flags);
3186    return result;
3187}
3188
3189template <>
3190uint32_t
3191fplibMax(uint32_t op1, uint32_t op2, FPSCR &fpscr)
3192{
3193    int mode = modeConv(fpscr);
3194    int flags = 0;
3195    int sgn1, exp1, sgn2, exp2;
3196    uint32_t mnt1, mnt2, x, result;
3197
3198    fp32_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
3199    fp32_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
3200
3201    if ((x = fp32_process_NaNs(op1, op2, mode, &flags))) {
3202        result = x;
3203    } else {
3204        result = ((sgn1 != sgn2 ? sgn2 : sgn1 ^ (op1 > op2)) ?
3205                  fp32_repack(sgn1, exp1, mnt1) :
3206                  fp32_repack(sgn2, exp2, mnt2));
3207    }
3208    set_fpscr0(fpscr, flags);
3209    return result;
3210}
3211
3212template <>
3213uint64_t
3214fplibMax(uint64_t op1, uint64_t op2, FPSCR &fpscr)
3215{
3216    int mode = modeConv(fpscr);
3217    int flags = 0;
3218    int sgn1, exp1, sgn2, exp2;
3219    uint64_t mnt1, mnt2, x, result;
3220
3221    fp64_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
3222    fp64_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
3223
3224    if ((x = fp64_process_NaNs(op1, op2, mode, &flags))) {
3225        result = x;
3226    } else {
3227        result = ((sgn1 != sgn2 ? sgn2 : sgn1 ^ (op1 > op2)) ?
3228                  fp64_repack(sgn1, exp1, mnt1) :
3229                  fp64_repack(sgn2, exp2, mnt2));
3230    }
3231    set_fpscr0(fpscr, flags);
3232    return result;
3233}
3234
3235template <>
3236uint16_t
3237fplibMaxNum(uint16_t op1, uint16_t op2, FPSCR &fpscr)
3238{
3239    fp16_minmaxnum(&op1, &op2, 1);
3240    return fplibMax<uint16_t>(op1, op2, fpscr);
3241}
3242
3243template <>
3244uint32_t
3245fplibMaxNum(uint32_t op1, uint32_t op2, FPSCR &fpscr)
3246{
3247    fp32_minmaxnum(&op1, &op2, 1);
3248    return fplibMax<uint32_t>(op1, op2, fpscr);
3249}
3250
3251template <>
3252uint64_t
3253fplibMaxNum(uint64_t op1, uint64_t op2, FPSCR &fpscr)
3254{
3255    fp64_minmaxnum(&op1, &op2, 1);
3256    return fplibMax<uint64_t>(op1, op2, fpscr);
3257}
3258
3259template <>
3260uint16_t
3261fplibMin(uint16_t op1, uint16_t op2, FPSCR &fpscr)
3262{
3263    int mode = modeConv(fpscr);
3264    int flags = 0;
3265    int sgn1, exp1, sgn2, exp2;
3266    uint16_t mnt1, mnt2, x, result;
3267
3268    fp16_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
3269    fp16_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
3270
3271    if ((x = fp16_process_NaNs(op1, op2, mode, &flags))) {
3272        result = x;
3273    } else {
3274        result = ((sgn1 != sgn2 ? sgn1 : sgn1 ^ (op1 < op2)) ?
3275                  fp16_repack(sgn1, exp1, mnt1) :
3276                  fp16_repack(sgn2, exp2, mnt2));
3277    }
3278    set_fpscr0(fpscr, flags);
3279    return result;
3280}
3281
3282template <>
3283uint32_t
3284fplibMin(uint32_t op1, uint32_t op2, FPSCR &fpscr)
3285{
3286    int mode = modeConv(fpscr);
3287    int flags = 0;
3288    int sgn1, exp1, sgn2, exp2;
3289    uint32_t mnt1, mnt2, x, result;
3290
3291    fp32_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
3292    fp32_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
3293
3294    if ((x = fp32_process_NaNs(op1, op2, mode, &flags))) {
3295        result = x;
3296    } else {
3297        result = ((sgn1 != sgn2 ? sgn1 : sgn1 ^ (op1 < op2)) ?
3298                  fp32_repack(sgn1, exp1, mnt1) :
3299                  fp32_repack(sgn2, exp2, mnt2));
3300    }
3301    set_fpscr0(fpscr, flags);
3302    return result;
3303}
3304
3305template <>
3306uint64_t
3307fplibMin(uint64_t op1, uint64_t op2, FPSCR &fpscr)
3308{
3309    int mode = modeConv(fpscr);
3310    int flags = 0;
3311    int sgn1, exp1, sgn2, exp2;
3312    uint64_t mnt1, mnt2, x, result;
3313
3314    fp64_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
3315    fp64_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
3316
3317    if ((x = fp64_process_NaNs(op1, op2, mode, &flags))) {
3318        result = x;
3319    } else {
3320        result = ((sgn1 != sgn2 ? sgn1 : sgn1 ^ (op1 < op2)) ?
3321                  fp64_repack(sgn1, exp1, mnt1) :
3322                  fp64_repack(sgn2, exp2, mnt2));
3323    }
3324    set_fpscr0(fpscr, flags);
3325    return result;
3326}
3327
3328template <>
3329uint16_t
3330fplibMinNum(uint16_t op1, uint16_t op2, FPSCR &fpscr)
3331{
3332    fp16_minmaxnum(&op1, &op2, 0);
3333    return fplibMin<uint16_t>(op1, op2, fpscr);
3334}
3335
3336template <>
3337uint32_t
3338fplibMinNum(uint32_t op1, uint32_t op2, FPSCR &fpscr)
3339{
3340    fp32_minmaxnum(&op1, &op2, 0);
3341    return fplibMin<uint32_t>(op1, op2, fpscr);
3342}
3343
3344template <>
3345uint64_t
3346fplibMinNum(uint64_t op1, uint64_t op2, FPSCR &fpscr)
3347{
3348    fp64_minmaxnum(&op1, &op2, 0);
3349    return fplibMin<uint64_t>(op1, op2, fpscr);
3350}
3351
3352template <>
3353uint16_t
3354fplibMul(uint16_t op1, uint16_t op2, FPSCR &fpscr)
3355{
3356    int flags = 0;
3357    uint16_t result = fp16_mul(op1, op2, modeConv(fpscr), &flags);
3358    set_fpscr0(fpscr, flags);
3359    return result;
3360}
3361
3362template <>
3363uint32_t
3364fplibMul(uint32_t op1, uint32_t op2, FPSCR &fpscr)
3365{
3366    int flags = 0;
3367    uint32_t result = fp32_mul(op1, op2, modeConv(fpscr), &flags);
3368    set_fpscr0(fpscr, flags);
3369    return result;
3370}
3371
3372template <>
3373uint64_t
3374fplibMul(uint64_t op1, uint64_t op2, FPSCR &fpscr)
3375{
3376    int flags = 0;
3377    uint64_t result = fp64_mul(op1, op2, modeConv(fpscr), &flags);
3378    set_fpscr0(fpscr, flags);
3379    return result;
3380}
3381
3382template <>
3383uint16_t
3384fplibMulX(uint16_t op1, uint16_t op2, FPSCR &fpscr)
3385{
3386    int mode = modeConv(fpscr);
3387    int flags = 0;
3388    int sgn1, exp1, sgn2, exp2;
3389    uint16_t mnt1, mnt2, result;
3390
3391    fp16_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
3392    fp16_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
3393
3394    result = fp16_process_NaNs(op1, op2, mode, &flags);
3395    if (!result) {
3396        if ((exp1 == FP16_EXP_INF && !mnt2) ||
3397            (exp2 == FP16_EXP_INF && !mnt1)) {
3398            result = fp16_FPTwo(sgn1 ^ sgn2);
3399        } else if (exp1 == FP16_EXP_INF || exp2 == FP16_EXP_INF) {
3400            result = fp16_infinity(sgn1 ^ sgn2);
3401        } else if (!mnt1 || !mnt2) {
3402            result = fp16_zero(sgn1 ^ sgn2);
3403        } else {
3404            result = fp16_mul(op1, op2, mode, &flags);
3405        }
3406    }
3407
3408    set_fpscr0(fpscr, flags);
3409
3410    return result;
3411}
3412
3413template <>
3414uint32_t
3415fplibMulX(uint32_t op1, uint32_t op2, FPSCR &fpscr)
3416{
3417    int mode = modeConv(fpscr);
3418    int flags = 0;
3419    int sgn1, exp1, sgn2, exp2;
3420    uint32_t mnt1, mnt2, result;
3421
3422    fp32_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
3423    fp32_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
3424
3425    result = fp32_process_NaNs(op1, op2, mode, &flags);
3426    if (!result) {
3427        if ((exp1 == FP32_EXP_INF && !mnt2) ||
3428            (exp2 == FP32_EXP_INF && !mnt1)) {
3429            result = fp32_FPTwo(sgn1 ^ sgn2);
3430        } else if (exp1 == FP32_EXP_INF || exp2 == FP32_EXP_INF) {
3431            result = fp32_infinity(sgn1 ^ sgn2);
3432        } else if (!mnt1 || !mnt2) {
3433            result = fp32_zero(sgn1 ^ sgn2);
3434        } else {
3435            result = fp32_mul(op1, op2, mode, &flags);
3436        }
3437    }
3438
3439    set_fpscr0(fpscr, flags);
3440
3441    return result;
3442}
3443
3444template <>
3445uint64_t
3446fplibMulX(uint64_t op1, uint64_t op2, FPSCR &fpscr)
3447{
3448    int mode = modeConv(fpscr);
3449    int flags = 0;
3450    int sgn1, exp1, sgn2, exp2;
3451    uint64_t mnt1, mnt2, result;
3452
3453    fp64_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
3454    fp64_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
3455
3456    result = fp64_process_NaNs(op1, op2, mode, &flags);
3457    if (!result) {
3458        if ((exp1 == FP64_EXP_INF && !mnt2) ||
3459            (exp2 == FP64_EXP_INF && !mnt1)) {
3460            result = fp64_FPTwo(sgn1 ^ sgn2);
3461        } else if (exp1 == FP64_EXP_INF || exp2 == FP64_EXP_INF) {
3462            result = fp64_infinity(sgn1 ^ sgn2);
3463        } else if (!mnt1 || !mnt2) {
3464            result = fp64_zero(sgn1 ^ sgn2);
3465        } else {
3466            result = fp64_mul(op1, op2, mode, &flags);
3467        }
3468    }
3469
3470    set_fpscr0(fpscr, flags);
3471
3472    return result;
3473}
3474
3475template <>
3476uint16_t
3477fplibNeg(uint16_t op)
3478{
3479    return op ^ 1ULL << (FP16_BITS - 1);
3480}
3481
3482template <>
3483uint32_t
3484fplibNeg(uint32_t op)
3485{
3486    return op ^ 1ULL << (FP32_BITS - 1);
3487}
3488
3489template <>
3490uint64_t
3491fplibNeg(uint64_t op)
3492{
3493    return op ^ 1ULL << (FP64_BITS - 1);
3494}
3495
3496static const uint8_t recip_sqrt_estimate[256] = {
3497    255, 253, 251, 249, 247, 245, 243, 242, 240, 238, 236, 234, 233, 231, 229, 228,
3498    226, 224, 223, 221, 219, 218, 216, 215, 213, 212, 210, 209, 207, 206, 204, 203,
3499    201, 200, 198, 197, 196, 194, 193, 192, 190, 189, 188, 186, 185, 184, 183, 181,
3500    180, 179, 178, 176, 175, 174, 173, 172, 170, 169, 168, 167, 166, 165, 164, 163,
3501    162, 160, 159, 158, 157, 156, 155, 154, 153, 152, 151, 150, 149, 148, 147, 146,
3502    145, 144, 143, 142, 141, 140, 140, 139, 138, 137, 136, 135, 134, 133, 132, 131,
3503    131, 130, 129, 128, 127, 126, 126, 125, 124, 123, 122, 121, 121, 120, 119, 118,
3504    118, 117, 116, 115, 114, 114, 113, 112, 111, 111, 110, 109, 109, 108, 107, 106,
3505    105, 104, 103, 101, 100,  99,  97,  96,  95,  93,  92,  91,  90,  88,  87,  86,
3506    85,  84,  82,  81,  80,  79,  78,  77,  76,  75,  74,  72,  71,  70,  69,  68,
3507    67,  66,  65,  64,  63,  62,  61,  60,  60,  59,  58,  57,  56,  55,  54,  53,
3508    52,  51,  51,  50,  49,  48,  47,  46,  46,  45,  44,  43,  42,  42,  41,  40,
3509    39,  38,  38,  37,  36,  35,  35,  34,  33,  33,  32,  31,  30,  30,  29,  28,
3510    28,  27,  26,  26,  25,  24,  24,  23,  22,  22,  21,  20,  20,  19,  19,  18,
3511    17,  17,  16,  16,  15,  14,  14,  13,  13,  12,  11,  11,  10,  10,   9,   9,
3512    8,   8,   7,   6,   6,   5,   5,   4,   4,   3,   3,   2,   2,   1,   1,   0
3513};
3514
3515template <>
3516uint16_t
3517fplibRSqrtEstimate(uint16_t op, FPSCR &fpscr)
3518{
3519    int mode = modeConv(fpscr);
3520    int flags = 0;
3521    int sgn, exp;
3522    uint16_t mnt, result;
3523
3524    fp16_unpack(&sgn, &exp, &mnt, op, mode, &flags);
3525
3526    if (fp16_is_NaN(exp, mnt)) {
3527        result = fp16_process_NaN(op, mode, &flags);
3528    } else if (!mnt) {
3529        result = fp16_infinity(sgn);
3530        flags |= FPLIB_DZC;
3531    } else if (sgn) {
3532        result = fp16_defaultNaN();
3533        flags |= FPLIB_IOC;
3534    } else if (exp == FP16_EXP_INF) {
3535        result = fp16_zero(0);
3536    } else {
3537        exp += FP16_EXP_BITS;
3538        mnt = fp16_normalise(mnt, &exp);
3539        mnt = recip_sqrt_estimate[(~exp & 1) << 7 |
3540                                  (mnt >> (FP16_BITS - 8) & 127)];
3541        result = fp16_pack(0, (3 * FP16_EXP_BIAS - exp - 1) >> 1,
3542                           mnt << (FP16_MANT_BITS - 8));
3543    }
3544
3545    set_fpscr0(fpscr, flags);
3546
3547    return result;
3548}
3549
3550template <>
3551uint32_t
3552fplibRSqrtEstimate(uint32_t op, FPSCR &fpscr)
3553{
3554    int mode = modeConv(fpscr);
3555    int flags = 0;
3556    int sgn, exp;
3557    uint32_t mnt, result;
3558
3559    fp32_unpack(&sgn, &exp, &mnt, op, mode, &flags);
3560
3561    if (fp32_is_NaN(exp, mnt)) {
3562        result = fp32_process_NaN(op, mode, &flags);
3563    } else if (!mnt) {
3564        result = fp32_infinity(sgn);
3565        flags |= FPLIB_DZC;
3566    } else if (sgn) {
3567        result = fp32_defaultNaN();
3568        flags |= FPLIB_IOC;
3569    } else if (exp == FP32_EXP_INF) {
3570        result = fp32_zero(0);
3571    } else {
3572        exp += FP32_EXP_BITS;
3573        mnt = fp32_normalise(mnt, &exp);
3574        mnt = recip_sqrt_estimate[(~exp & 1) << 7 |
3575                                  (mnt >> (FP32_BITS - 8) & 127)];
3576        result = fp32_pack(0, (3 * FP32_EXP_BIAS - exp - 1) >> 1,
3577                           mnt << (FP32_MANT_BITS - 8));
3578    }
3579
3580    set_fpscr0(fpscr, flags);
3581
3582    return result;
3583}
3584
3585template <>
3586uint64_t
3587fplibRSqrtEstimate(uint64_t op, FPSCR &fpscr)
3588{
3589    int mode = modeConv(fpscr);
3590    int flags = 0;
3591    int sgn, exp;
3592    uint64_t mnt, result;
3593
3594    fp64_unpack(&sgn, &exp, &mnt, op, mode, &flags);
3595
3596    if (fp64_is_NaN(exp, mnt)) {
3597        result = fp64_process_NaN(op, mode, &flags);
3598    } else if (!mnt) {
3599        result = fp64_infinity(sgn);
3600        flags |= FPLIB_DZC;
3601    } else if (sgn) {
3602        result = fp64_defaultNaN();
3603        flags |= FPLIB_IOC;
3604    } else if (exp == FP64_EXP_INF) {
3605        result = fp32_zero(0);
3606    } else {
3607        exp += FP64_EXP_BITS;
3608        mnt = fp64_normalise(mnt, &exp);
3609        mnt = recip_sqrt_estimate[(~exp & 1) << 7 |
3610                                  (mnt >> (FP64_BITS - 8) & 127)];
3611        result = fp64_pack(0, (3 * FP64_EXP_BIAS - exp - 1) >> 1,
3612                           mnt << (FP64_MANT_BITS - 8));
3613    }
3614
3615    set_fpscr0(fpscr, flags);
3616
3617    return result;
3618}
3619
3620template <>
3621uint16_t
3622fplibRSqrtStepFused(uint16_t op1, uint16_t op2, FPSCR &fpscr)
3623{
3624    int mode = modeConv(fpscr);
3625    int flags = 0;
3626    int sgn1, exp1, sgn2, exp2;
3627    uint16_t mnt1, mnt2, result;
3628
3629    op1 = fplibNeg<uint16_t>(op1);
3630    fp16_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
3631    fp16_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
3632
3633    result = fp16_process_NaNs(op1, op2, mode, &flags);
3634    if (!result) {
3635        if ((exp1 == FP16_EXP_INF && !mnt2) ||
3636            (exp2 == FP16_EXP_INF && !mnt1)) {
3637            result = fp16_FPOnePointFive(0);
3638        } else if (exp1 == FP16_EXP_INF || exp2 == FP16_EXP_INF) {
3639            result = fp16_infinity(sgn1 ^ sgn2);
3640        } else {
3641            result = fp16_muladd(fp16_FPThree(0), op1, op2, -1, mode, &flags);
3642        }
3643    }
3644
3645    set_fpscr0(fpscr, flags);
3646
3647    return result;
3648}
3649
3650template <>
3651uint32_t
3652fplibRSqrtStepFused(uint32_t op1, uint32_t op2, FPSCR &fpscr)
3653{
3654    int mode = modeConv(fpscr);
3655    int flags = 0;
3656    int sgn1, exp1, sgn2, exp2;
3657    uint32_t mnt1, mnt2, result;
3658
3659    op1 = fplibNeg<uint32_t>(op1);
3660    fp32_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
3661    fp32_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
3662
3663    result = fp32_process_NaNs(op1, op2, mode, &flags);
3664    if (!result) {
3665        if ((exp1 == FP32_EXP_INF && !mnt2) ||
3666            (exp2 == FP32_EXP_INF && !mnt1)) {
3667            result = fp32_FPOnePointFive(0);
3668        } else if (exp1 == FP32_EXP_INF || exp2 == FP32_EXP_INF) {
3669            result = fp32_infinity(sgn1 ^ sgn2);
3670        } else {
3671            result = fp32_muladd(fp32_FPThree(0), op1, op2, -1, mode, &flags);
3672        }
3673    }
3674
3675    set_fpscr0(fpscr, flags);
3676
3677    return result;
3678}
3679
3680template <>
3681uint64_t
3682fplibRSqrtStepFused(uint64_t op1, uint64_t op2, FPSCR &fpscr)
3683{
3684    int mode = modeConv(fpscr);
3685    int flags = 0;
3686    int sgn1, exp1, sgn2, exp2;
3687    uint64_t mnt1, mnt2, result;
3688
3689    op1 = fplibNeg<uint64_t>(op1);
3690    fp64_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
3691    fp64_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
3692
3693    result = fp64_process_NaNs(op1, op2, mode, &flags);
3694    if (!result) {
3695        if ((exp1 == FP64_EXP_INF && !mnt2) ||
3696            (exp2 == FP64_EXP_INF && !mnt1)) {
3697            result = fp64_FPOnePointFive(0);
3698        } else if (exp1 == FP64_EXP_INF || exp2 == FP64_EXP_INF) {
3699            result = fp64_infinity(sgn1 ^ sgn2);
3700        } else {
3701            result = fp64_muladd(fp64_FPThree(0), op1, op2, -1, mode, &flags);
3702        }
3703    }
3704
3705    set_fpscr0(fpscr, flags);
3706
3707    return result;
3708}
3709
3710template <>
3711uint16_t
3712fplibRecipEstimate(uint16_t op, FPSCR &fpscr)
3713{
3714    int mode = modeConv(fpscr);
3715    int flags = 0;
3716    int sgn, exp;
3717    uint16_t mnt, result;
3718
3719    fp16_unpack(&sgn, &exp, &mnt, op, mode, &flags);
3720
3721    if (fp16_is_NaN(exp, mnt)) {
3722        result = fp16_process_NaN(op, mode, &flags);
3723    } else if (exp == FP16_EXP_INF) {
3724        result = fp16_zero(sgn);
3725    } else if (!mnt) {
3726        result = fp16_infinity(sgn);
3727        flags |= FPLIB_DZC;
3728    } else if (!((uint16_t)(op << 1) >> (FP16_MANT_BITS - 1))) {
3729        bool overflow_to_inf = false;
3730        switch (FPCRRounding(fpscr)) {
3731          case FPRounding_TIEEVEN:
3732            overflow_to_inf = true;
3733            break;
3734          case FPRounding_POSINF:
3735            overflow_to_inf = !sgn;
3736            break;
3737          case FPRounding_NEGINF:
3738            overflow_to_inf = sgn;
3739            break;
3740          case FPRounding_ZERO:
3741            overflow_to_inf = false;
3742            break;
3743          default:
3744            panic("Unrecognized FP rounding mode");
3745        }
3746        result = overflow_to_inf ? fp16_infinity(sgn) : fp16_max_normal(sgn);
3747        flags |= FPLIB_OFC | FPLIB_IXC;
3748    } else if (fpscr.fz16 && exp >= 2 * FP16_EXP_BIAS - 1) {
3749        result = fp16_zero(sgn);
3750        flags |= FPLIB_UFC;
3751    } else {
3752        exp += FP16_EXP_BITS;
3753        mnt = fp16_normalise(mnt, &exp);
3754        int result_exp = 2 * FP16_EXP_BIAS - 1 - exp;
3755        uint16_t fraction = (((uint32_t)1 << 19) /
3756                             (mnt >> (FP16_BITS - 10) | 1) + 1) >> 1;
3757        fraction <<= FP16_MANT_BITS - 8;
3758        if (result_exp == 0) {
3759            fraction >>= 1;
3760        } else if (result_exp == -1) {
3761            fraction >>= 2;
3762            result_exp = 0;
3763        }
3764        result = fp16_pack(sgn, result_exp, fraction);
3765    }
3766
3767    set_fpscr0(fpscr, flags);
3768
3769    return result;
3770}
3771
3772template <>
3773uint32_t
3774fplibRecipEstimate(uint32_t op, FPSCR &fpscr)
3775{
3776    int mode = modeConv(fpscr);
3777    int flags = 0;
3778    int sgn, exp;
3779    uint32_t mnt, result;
3780
3781    fp32_unpack(&sgn, &exp, &mnt, op, mode, &flags);
3782
3783    if (fp32_is_NaN(exp, mnt)) {
3784        result = fp32_process_NaN(op, mode, &flags);
3785    } else if (exp == FP32_EXP_INF) {
3786        result = fp32_zero(sgn);
3787    } else if (!mnt) {
3788        result = fp32_infinity(sgn);
3789        flags |= FPLIB_DZC;
3790    } else if (!((uint32_t)(op << 1) >> (FP32_MANT_BITS - 1))) {
3791        bool overflow_to_inf = false;
3792        switch (FPCRRounding(fpscr)) {
3793          case FPRounding_TIEEVEN:
3794            overflow_to_inf = true;
3795            break;
3796          case FPRounding_POSINF:
3797            overflow_to_inf = !sgn;
3798            break;
3799          case FPRounding_NEGINF:
3800            overflow_to_inf = sgn;
3801            break;
3802          case FPRounding_ZERO:
3803            overflow_to_inf = false;
3804            break;
3805          default:
3806            panic("Unrecognized FP rounding mode");
3807        }
3808        result = overflow_to_inf ? fp32_infinity(sgn) : fp32_max_normal(sgn);
3809        flags |= FPLIB_OFC | FPLIB_IXC;
3810    } else if (fpscr.fz && exp >= 2 * FP32_EXP_BIAS - 1) {
3811        result = fp32_zero(sgn);
3812        flags |= FPLIB_UFC;
3813    } else {
3814        exp += FP32_EXP_BITS;
3815        mnt = fp32_normalise(mnt, &exp);
3816        int result_exp = 2 * FP32_EXP_BIAS - 1 - exp;
3817        uint32_t fraction = (((uint32_t)1 << 19) /
3818                             (mnt >> (FP32_BITS - 10) | 1) + 1) >> 1;
3819        fraction <<= FP32_MANT_BITS - 8;
3820        if (result_exp == 0) {
3821            fraction >>= 1;
3822        } else if (result_exp == -1) {
3823            fraction >>= 2;
3824            result_exp = 0;
3825        }
3826        result = fp32_pack(sgn, result_exp, fraction);
3827    }
3828
3829    set_fpscr0(fpscr, flags);
3830
3831    return result;
3832}
3833
3834template <>
3835uint64_t
3836fplibRecipEstimate(uint64_t op, FPSCR &fpscr)
3837{
3838    int mode = modeConv(fpscr);
3839    int flags = 0;
3840    int sgn, exp;
3841    uint64_t mnt, result;
3842
3843    fp64_unpack(&sgn, &exp, &mnt, op, mode, &flags);
3844
3845    if (fp64_is_NaN(exp, mnt)) {
3846        result = fp64_process_NaN(op, mode, &flags);
3847    } else if (exp == FP64_EXP_INF) {
3848        result = fp64_zero(sgn);
3849    } else if (!mnt) {
3850        result = fp64_infinity(sgn);
3851        flags |= FPLIB_DZC;
3852    } else if (!((uint64_t)(op << 1) >> (FP64_MANT_BITS - 1))) {
3853        bool overflow_to_inf = false;
3854        switch (FPCRRounding(fpscr)) {
3855          case FPRounding_TIEEVEN:
3856            overflow_to_inf = true;
3857            break;
3858          case FPRounding_POSINF:
3859            overflow_to_inf = !sgn;
3860            break;
3861          case FPRounding_NEGINF:
3862            overflow_to_inf = sgn;
3863            break;
3864          case FPRounding_ZERO:
3865            overflow_to_inf = false;
3866            break;
3867          default:
3868            panic("Unrecognized FP rounding mode");
3869        }
3870        result = overflow_to_inf ? fp64_infinity(sgn) : fp64_max_normal(sgn);
3871        flags |= FPLIB_OFC | FPLIB_IXC;
3872    } else if (fpscr.fz && exp >= 2 * FP64_EXP_BIAS - 1) {
3873        result = fp64_zero(sgn);
3874        flags |= FPLIB_UFC;
3875    } else {
3876        exp += FP64_EXP_BITS;
3877        mnt = fp64_normalise(mnt, &exp);
3878        int result_exp = 2 * FP64_EXP_BIAS - 1 - exp;
3879        uint64_t fraction = (((uint32_t)1 << 19) /
3880                             (mnt >> (FP64_BITS - 10) | 1) + 1) >> 1;
3881        fraction <<= FP64_MANT_BITS - 8;
3882        if (result_exp == 0) {
3883            fraction >>= 1;
3884        } else if (result_exp == -1) {
3885            fraction >>= 2;
3886            result_exp = 0;
3887        }
3888        result = fp64_pack(sgn, result_exp, fraction);
3889    }
3890
3891    set_fpscr0(fpscr, flags);
3892
3893    return result;
3894}
3895
3896template <>
3897uint16_t
3898fplibRecipStepFused(uint16_t op1, uint16_t op2, FPSCR &fpscr)
3899{
3900    int mode = modeConv(fpscr);
3901    int flags = 0;
3902    int sgn1, exp1, sgn2, exp2;
3903    uint16_t mnt1, mnt2, result;
3904
3905    op1 = fplibNeg<uint16_t>(op1);
3906    fp16_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
3907    fp16_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
3908
3909    result = fp16_process_NaNs(op1, op2, mode, &flags);
3910    if (!result) {
3911        if ((exp1 == FP16_EXP_INF && !mnt2) ||
3912            (exp2 == FP16_EXP_INF && !mnt1)) {
3913            result = fp16_FPTwo(0);
3914        } else if (exp1 == FP16_EXP_INF || exp2 == FP16_EXP_INF) {
3915            result = fp16_infinity(sgn1 ^ sgn2);
3916        } else {
3917            result = fp16_muladd(fp16_FPTwo(0), op1, op2, 0, mode, &flags);
3918        }
3919    }
3920
3921    set_fpscr0(fpscr, flags);
3922
3923    return result;
3924}
3925
3926template <>
3927uint32_t
3928fplibRecipStepFused(uint32_t op1, uint32_t op2, FPSCR &fpscr)
3929{
3930    int mode = modeConv(fpscr);
3931    int flags = 0;
3932    int sgn1, exp1, sgn2, exp2;
3933    uint32_t mnt1, mnt2, result;
3934
3935    op1 = fplibNeg<uint32_t>(op1);
3936    fp32_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
3937    fp32_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
3938
3939    result = fp32_process_NaNs(op1, op2, mode, &flags);
3940    if (!result) {
3941        if ((exp1 == FP32_EXP_INF && !mnt2) ||
3942            (exp2 == FP32_EXP_INF && !mnt1)) {
3943            result = fp32_FPTwo(0);
3944        } else if (exp1 == FP32_EXP_INF || exp2 == FP32_EXP_INF) {
3945            result = fp32_infinity(sgn1 ^ sgn2);
3946        } else {
3947            result = fp32_muladd(fp32_FPTwo(0), op1, op2, 0, mode, &flags);
3948        }
3949    }
3950
3951    set_fpscr0(fpscr, flags);
3952
3953    return result;
3954}
3955
3956template <>
3957uint64_t
3958fplibRecipStepFused(uint64_t op1, uint64_t op2, FPSCR &fpscr)
3959{
3960    int mode = modeConv(fpscr);
3961    int flags = 0;
3962    int sgn1, exp1, sgn2, exp2;
3963    uint64_t mnt1, mnt2, result;
3964
3965    op1 = fplibNeg<uint64_t>(op1);
3966    fp64_unpack(&sgn1, &exp1, &mnt1, op1, mode, &flags);
3967    fp64_unpack(&sgn2, &exp2, &mnt2, op2, mode, &flags);
3968
3969    result = fp64_process_NaNs(op1, op2, mode, &flags);
3970    if (!result) {
3971        if ((exp1 == FP64_EXP_INF && !mnt2) ||
3972            (exp2 == FP64_EXP_INF && !mnt1)) {
3973            result = fp64_FPTwo(0);
3974        } else if (exp1 == FP64_EXP_INF || exp2 == FP64_EXP_INF) {
3975            result = fp64_infinity(sgn1 ^ sgn2);
3976        } else {
3977            result = fp64_muladd(fp64_FPTwo(0), op1, op2, 0, mode, &flags);
3978        }
3979    }
3980
3981    set_fpscr0(fpscr, flags);
3982
3983    return result;
3984}
3985
3986template <>
3987uint16_t
3988fplibRecpX(uint16_t op, FPSCR &fpscr)
3989{
3990    int mode = modeConv(fpscr);
3991    int flags = 0;
3992    int sgn, exp;
3993    uint16_t mnt, result;
3994
3995    fp16_unpack(&sgn, &exp, &mnt, op, mode, &flags);
3996
3997    if (fp16_is_NaN(exp, mnt)) {
3998        result = fp16_process_NaN(op, mode, &flags);
3999    }
4000    else {
4001        if (!mnt) { // Zero and denormals
4002            result = fp16_pack(sgn, FP16_EXP_INF - 1, 0);
4003        } else { // Infinities and normals
4004            result = fp16_pack(sgn, exp ^ FP16_EXP_INF, 0);
4005        }
4006    }
4007
4008    set_fpscr0(fpscr, flags);
4009
4010    return result;
4011}
4012
4013template <>
4014uint32_t
4015fplibRecpX(uint32_t op, FPSCR &fpscr)
4016{
4017    int mode = modeConv(fpscr);
4018    int flags = 0;
4019    int sgn, exp;
4020    uint32_t mnt, result;
4021
4022    fp32_unpack(&sgn, &exp, &mnt, op, mode, &flags);
4023
4024    if (fp32_is_NaN(exp, mnt)) {
4025        result = fp32_process_NaN(op, mode, &flags);
4026    }
4027    else {
4028        if (!mnt) { // Zero and denormals
4029            result = fp32_pack(sgn, FP32_EXP_INF - 1, 0);
4030        } else { // Infinities and normals
4031            result = fp32_pack(sgn, exp ^ FP32_EXP_INF, 0);
4032        }
4033    }
4034
4035    set_fpscr0(fpscr, flags);
4036
4037    return result;
4038}
4039
4040template <>
4041uint64_t
4042fplibRecpX(uint64_t op, FPSCR &fpscr)
4043{
4044    int mode = modeConv(fpscr);
4045    int flags = 0;
4046    int sgn, exp;
4047    uint64_t mnt, result;
4048
4049    fp64_unpack(&sgn, &exp, &mnt, op, mode, &flags);
4050
4051    if (fp64_is_NaN(exp, mnt)) {
4052        result = fp64_process_NaN(op, mode, &flags);
4053    }
4054    else {
4055        if (!mnt) { // Zero and denormals
4056            result = fp64_pack(sgn, FP64_EXP_INF - 1, 0);
4057        } else { // Infinities and normals
4058            result = fp64_pack(sgn, exp ^ FP64_EXP_INF, 0);
4059        }
4060    }
4061
4062    set_fpscr0(fpscr, flags);
4063
4064    return result;
4065}
4066
4067template <>
4068uint16_t
4069fplibRoundInt(uint16_t op, FPRounding rounding, bool exact, FPSCR &fpscr)
4070{
4071    int expint = FP16_EXP_BIAS + FP16_MANT_BITS;
4072    int mode = modeConv(fpscr);
4073    int flags = 0;
4074    int sgn, exp;
4075    uint16_t mnt, result;
4076
4077    // Unpack using FPCR to determine if subnormals are flushed-to-zero:
4078    fp16_unpack(&sgn, &exp, &mnt, op, mode, &flags);
4079
4080    // Handle NaNs, infinities and zeroes:
4081    if (fp16_is_NaN(exp, mnt)) {
4082        result = fp16_process_NaN(op, mode, &flags);
4083    } else if (exp == FP16_EXP_INF) {
4084        result = fp16_infinity(sgn);
4085    } else if (!mnt) {
4086        result = fp16_zero(sgn);
4087    } else if (exp >= expint) {
4088        // There are no fractional bits
4089        result = op;
4090    } else {
4091        // Truncate towards zero:
4092        uint16_t x = expint - exp >= FP16_BITS ? 0 : mnt >> (expint - exp);
4093        int err = exp < expint - FP16_BITS ? 1 :
4094            ((mnt << 1 >> (expint - exp - 1) & 3) |
4095             ((uint16_t)(mnt << 2 << (FP16_BITS + exp - expint)) != 0));
4096        switch (rounding) {
4097          case FPRounding_TIEEVEN:
4098            x += (err == 3 || (err == 2 && (x & 1)));
4099            break;
4100          case FPRounding_POSINF:
4101            x += err && !sgn;
4102            break;
4103          case FPRounding_NEGINF:
4104            x += err && sgn;
4105            break;
4106          case FPRounding_ZERO:
4107            break;
4108          case FPRounding_TIEAWAY:
4109            x += err >> 1;
4110            break;
4111          default:
4112            panic("Unrecognized FP rounding mode");
4113        }
4114
4115        if (x == 0) {
4116            result = fp16_zero(sgn);
4117        } else {
4118            exp = expint;
4119            mnt = fp16_normalise(x, &exp);
4120            result = fp16_pack(sgn, exp + FP16_EXP_BITS, mnt >> FP16_EXP_BITS);
4121        }
4122
4123        if (err && exact)
4124            flags |= FPLIB_IXC;
4125    }
4126
4127    set_fpscr0(fpscr, flags);
4128
4129    return result;
4130}
4131
4132template <>
4133uint32_t
4134fplibRoundInt(uint32_t op, FPRounding rounding, bool exact, FPSCR &fpscr)
4135{
4136    int expint = FP32_EXP_BIAS + FP32_MANT_BITS;
4137    int mode = modeConv(fpscr);
4138    int flags = 0;
4139    int sgn, exp;
4140    uint32_t mnt, result;
4141
4142    // Unpack using FPCR to determine if subnormals are flushed-to-zero:
4143    fp32_unpack(&sgn, &exp, &mnt, op, mode, &flags);
4144
4145    // Handle NaNs, infinities and zeroes:
4146    if (fp32_is_NaN(exp, mnt)) {
4147        result = fp32_process_NaN(op, mode, &flags);
4148    } else if (exp == FP32_EXP_INF) {
4149        result = fp32_infinity(sgn);
4150    } else if (!mnt) {
4151        result = fp32_zero(sgn);
4152    } else if (exp >= expint) {
4153        // There are no fractional bits
4154        result = op;
4155    } else {
4156        // Truncate towards zero:
4157        uint32_t x = expint - exp >= FP32_BITS ? 0 : mnt >> (expint - exp);
4158        int err = exp < expint - FP32_BITS ? 1 :
4159            ((mnt << 1 >> (expint - exp - 1) & 3) |
4160             ((uint32_t)(mnt << 2 << (FP32_BITS + exp - expint)) != 0));
4161        switch (rounding) {
4162          case FPRounding_TIEEVEN:
4163            x += (err == 3 || (err == 2 && (x & 1)));
4164            break;
4165          case FPRounding_POSINF:
4166            x += err && !sgn;
4167            break;
4168          case FPRounding_NEGINF:
4169            x += err && sgn;
4170            break;
4171          case FPRounding_ZERO:
4172            break;
4173          case FPRounding_TIEAWAY:
4174            x += err >> 1;
4175            break;
4176          default:
4177            panic("Unrecognized FP rounding mode");
4178        }
4179
4180        if (x == 0) {
4181            result = fp32_zero(sgn);
4182        } else {
4183            exp = expint;
4184            mnt = fp32_normalise(x, &exp);
4185            result = fp32_pack(sgn, exp + FP32_EXP_BITS, mnt >> FP32_EXP_BITS);
4186        }
4187
4188        if (err && exact)
4189            flags |= FPLIB_IXC;
4190    }
4191
4192    set_fpscr0(fpscr, flags);
4193
4194    return result;
4195}
4196
4197template <>
4198uint64_t
4199fplibRoundInt(uint64_t op, FPRounding rounding, bool exact, FPSCR &fpscr)
4200{
4201    int expint = FP64_EXP_BIAS + FP64_MANT_BITS;
4202    int mode = modeConv(fpscr);
4203    int flags = 0;
4204    int sgn, exp;
4205    uint64_t mnt, result;
4206
4207    // Unpack using FPCR to determine if subnormals are flushed-to-zero:
4208    fp64_unpack(&sgn, &exp, &mnt, op, mode, &flags);
4209
4210    // Handle NaNs, infinities and zeroes:
4211    if (fp64_is_NaN(exp, mnt)) {
4212        result = fp64_process_NaN(op, mode, &flags);
4213    } else if (exp == FP64_EXP_INF) {
4214        result = fp64_infinity(sgn);
4215    } else if (!mnt) {
4216        result = fp64_zero(sgn);
4217    } else if (exp >= expint) {
4218        // There are no fractional bits
4219        result = op;
4220    } else {
4221        // Truncate towards zero:
4222        uint64_t x = expint - exp >= FP64_BITS ? 0 : mnt >> (expint - exp);
4223        int err = exp < expint - FP64_BITS ? 1 :
4224            ((mnt << 1 >> (expint - exp - 1) & 3) |
4225             ((uint64_t)(mnt << 2 << (FP64_BITS + exp - expint)) != 0));
4226        switch (rounding) {
4227          case FPRounding_TIEEVEN:
4228            x += (err == 3 || (err == 2 && (x & 1)));
4229            break;
4230          case FPRounding_POSINF:
4231            x += err && !sgn;
4232            break;
4233          case FPRounding_NEGINF:
4234            x += err && sgn;
4235            break;
4236          case FPRounding_ZERO:
4237            break;
4238          case FPRounding_TIEAWAY:
4239            x += err >> 1;
4240            break;
4241          default:
4242            panic("Unrecognized FP rounding mode");
4243        }
4244
4245        if (x == 0) {
4246            result = fp64_zero(sgn);
4247        } else {
4248            exp = expint;
4249            mnt = fp64_normalise(x, &exp);
4250            result = fp64_pack(sgn, exp + FP64_EXP_BITS, mnt >> FP64_EXP_BITS);
4251        }
4252
4253        if (err && exact)
4254            flags |= FPLIB_IXC;
4255    }
4256
4257    set_fpscr0(fpscr, flags);
4258
4259    return result;
4260}
4261
4262template <>
4263uint16_t
4264fplibScale(uint16_t op1, uint16_t op2, FPSCR &fpscr)
4265{
4266    int flags = 0;
4267    uint16_t result = fp16_scale(op1, (int16_t)op2, modeConv(fpscr), &flags);
4268    set_fpscr0(fpscr, flags);
4269    return result;
4270}
4271
4272template <>
4273uint32_t
4274fplibScale(uint32_t op1, uint32_t op2, FPSCR &fpscr)
4275{
4276    int flags = 0;
4277    uint32_t result = fp32_scale(op1, (int32_t)op2, modeConv(fpscr), &flags);
4278    set_fpscr0(fpscr, flags);
4279    return result;
4280}
4281
4282template <>
4283uint64_t
4284fplibScale(uint64_t op1, uint64_t op2, FPSCR &fpscr)
4285{
4286    int flags = 0;
4287    uint64_t result = fp64_scale(op1, (int64_t)op2, modeConv(fpscr), &flags);
4288    set_fpscr0(fpscr, flags);
4289    return result;
4290}
4291
4292template <>
4293uint16_t
4294fplibSqrt(uint16_t op, FPSCR &fpscr)
4295{
4296    int flags = 0;
4297    uint16_t result = fp16_sqrt(op, modeConv(fpscr), &flags);
4298    set_fpscr0(fpscr, flags);
4299    return result;
4300}
4301
4302template <>
4303uint32_t
4304fplibSqrt(uint32_t op, FPSCR &fpscr)
4305{
4306    int flags = 0;
4307    uint32_t result = fp32_sqrt(op, modeConv(fpscr), &flags);
4308    set_fpscr0(fpscr, flags);
4309    return result;
4310}
4311
4312template <>
4313uint64_t
4314fplibSqrt(uint64_t op, FPSCR &fpscr)
4315{
4316    int flags = 0;
4317    uint64_t result = fp64_sqrt(op, modeConv(fpscr), &flags);
4318    set_fpscr0(fpscr, flags);
4319    return result;
4320}
4321
4322template <>
4323uint16_t
4324fplibSub(uint16_t op1, uint16_t op2, FPSCR &fpscr)
4325{
4326    int flags = 0;
4327    uint16_t result = fp16_add(op1, op2, 1, modeConv(fpscr), &flags);
4328    set_fpscr0(fpscr, flags);
4329    return result;
4330}
4331
4332template <>
4333uint32_t
4334fplibSub(uint32_t op1, uint32_t op2, FPSCR &fpscr)
4335{
4336    int flags = 0;
4337    uint32_t result = fp32_add(op1, op2, 1, modeConv(fpscr), &flags);
4338    set_fpscr0(fpscr, flags);
4339    return result;
4340}
4341
4342template <>
4343uint64_t
4344fplibSub(uint64_t op1, uint64_t op2, FPSCR &fpscr)
4345{
4346    int flags = 0;
4347    uint64_t result = fp64_add(op1, op2, 1, modeConv(fpscr), &flags);
4348    set_fpscr0(fpscr, flags);
4349    return result;
4350}
4351
4352template <>
4353uint16_t
4354fplibTrigMulAdd(uint8_t coeff_index, uint16_t op1, uint16_t op2, FPSCR &fpscr)
4355{
4356    static uint16_t coeff[2][8] = {
4357        {
4358            0x3c00,
4359            0xb155,
4360            0x2030,
4361            0x0000,
4362            0x0000,
4363            0x0000,
4364            0x0000,
4365            0x0000,
4366        },
4367        {
4368            0x3c00,
4369            0xb800,
4370            0x293a,
4371            0x0000,
4372            0x0000,
4373            0x0000,
4374            0x0000,
4375            0x0000
4376        }
4377    };
4378    int flags = 0;
4379    uint16_t result =
4380        fp16_muladd(coeff[op2 >> (FP16_BITS - 1)][coeff_index], op1,
4381                    fplibAbs(op2), 0, modeConv(fpscr), &flags);
4382    set_fpscr0(fpscr, flags);
4383    return result;
4384}
4385
4386template <>
4387uint32_t
4388fplibTrigMulAdd(uint8_t coeff_index, uint32_t op1, uint32_t op2, FPSCR &fpscr)
4389{
4390    static uint32_t coeff[2][8] = {
4391        {
4392            0x3f800000,
4393            0xbe2aaaab,
4394            0x3c088886,
4395            0xb95008b9,
4396            0x36369d6d,
4397            0x00000000,
4398            0x00000000,
4399            0x00000000
4400        },
4401        {
4402            0x3f800000,
4403            0xbf000000,
4404            0x3d2aaaa6,
4405            0xbab60705,
4406            0x37cd37cc,
4407            0x00000000,
4408            0x00000000,
4409            0x00000000
4410        }
4411    };
4412    int flags = 0;
4413    uint32_t result =
4414        fp32_muladd(coeff[op2 >> (FP32_BITS - 1)][coeff_index], op1,
4415                    fplibAbs(op2), 0, modeConv(fpscr), &flags);
4416    set_fpscr0(fpscr, flags);
4417    return result;
4418}
4419
4420template <>
4421uint64_t
4422fplibTrigMulAdd(uint8_t coeff_index, uint64_t op1, uint64_t op2, FPSCR &fpscr)
4423{
4424    static uint64_t coeff[2][8] = {
4425        {
4426            0x3ff0000000000000ULL,
4427            0xbfc5555555555543ULL,
4428            0x3f8111111110f30cULL,
4429            0xbf2a01a019b92fc6ULL,
4430            0x3ec71de351f3d22bULL,
4431            0xbe5ae5e2b60f7b91ULL,
4432            0x3de5d8408868552fULL,
4433            0x0000000000000000ULL
4434        },
4435        {
4436            0x3ff0000000000000ULL,
4437            0xbfe0000000000000ULL,
4438            0x3fa5555555555536ULL,
4439            0xbf56c16c16c13a0bULL,
4440            0x3efa01a019b1e8d8ULL,
4441            0xbe927e4f7282f468ULL,
4442            0x3e21ee96d2641b13ULL,
4443            0xbda8f76380fbb401ULL
4444        }
4445    };
4446    int flags = 0;
4447    uint64_t result =
4448        fp64_muladd(coeff[op2 >> (FP64_BITS - 1)][coeff_index], op1,
4449                    fplibAbs(op2), 0, modeConv(fpscr), &flags);
4450    set_fpscr0(fpscr, flags);
4451    return result;
4452}
4453
4454template <>
4455uint16_t
4456fplibTrigSMul(uint16_t op1, uint16_t op2, FPSCR &fpscr)
4457{
4458    int flags = 0;
4459    int sgn, exp;
4460    uint16_t mnt;
4461
4462    int mode = modeConv(fpscr);
4463    uint16_t result = fp16_mul(op1, op1, mode, &flags);
4464    set_fpscr0(fpscr, flags);
4465
4466    fp16_unpack(&sgn, &exp, &mnt, result, mode, &flags);
4467    if (!fp16_is_NaN(exp, mnt)) {
4468        result = (result & ~(1ULL << (FP16_BITS - 1))) |
4469            op2 << (FP16_BITS - 1);
4470    }
4471    return result;
4472}
4473
4474template <>
4475uint32_t
4476fplibTrigSMul(uint32_t op1, uint32_t op2, FPSCR &fpscr)
4477{
4478    int flags = 0;
4479    int sgn, exp;
4480    uint32_t mnt;
4481
4482    int mode = modeConv(fpscr);
4483    uint32_t result = fp32_mul(op1, op1, mode, &flags);
4484    set_fpscr0(fpscr, flags);
4485
4486    fp32_unpack(&sgn, &exp, &mnt, result, mode, &flags);
4487    if (!fp32_is_NaN(exp, mnt)) {
4488        result = (result & ~(1ULL << (FP32_BITS - 1))) | op2 << (FP32_BITS - 1);
4489    }
4490    return result;
4491}
4492
4493template <>
4494uint64_t
4495fplibTrigSMul(uint64_t op1, uint64_t op2, FPSCR &fpscr)
4496{
4497    int flags = 0;
4498    int sgn, exp;
4499    uint64_t mnt;
4500
4501    int mode = modeConv(fpscr);
4502    uint64_t result = fp64_mul(op1, op1, mode, &flags);
4503    set_fpscr0(fpscr, flags);
4504
4505    fp64_unpack(&sgn, &exp, &mnt, result, mode, &flags);
4506    if (!fp64_is_NaN(exp, mnt)) {
4507        result = (result & ~(1ULL << (FP64_BITS - 1))) | op2 << (FP64_BITS - 1);
4508    }
4509    return result;
4510}
4511
4512template <>
4513uint16_t
4514fplibTrigSSel(uint16_t op1, uint16_t op2, FPSCR &fpscr)
4515{
4516    static constexpr uint16_t fpOne =
4517        (uint16_t)FP16_EXP_BIAS << FP16_MANT_BITS; // 1.0
4518    if (op2 & 1)
4519        op1 = fpOne;
4520    return op1 ^ ((op2 >> 1) << (FP16_BITS - 1));
4521}
4522
4523template <>
4524uint32_t
4525fplibTrigSSel(uint32_t op1, uint32_t op2, FPSCR &fpscr)
4526{
4527    static constexpr uint32_t fpOne =
4528        (uint32_t)FP32_EXP_BIAS << FP32_MANT_BITS; // 1.0
4529    if (op2 & 1)
4530        op1 = fpOne;
4531    return op1 ^ ((op2 >> 1) << (FP32_BITS - 1));
4532}
4533
4534template <>
4535uint64_t
4536fplibTrigSSel(uint64_t op1, uint64_t op2, FPSCR &fpscr)
4537{
4538    static constexpr uint64_t fpOne =
4539        (uint64_t)FP64_EXP_BIAS << FP64_MANT_BITS; // 1.0
4540    if (op2 & 1)
4541        op1 = fpOne;
4542    return op1 ^ ((op2 >> 1) << (FP64_BITS - 1));
4543}
4544
4545static uint64_t
4546FPToFixed_64(int sgn, int exp, uint64_t mnt, bool u, FPRounding rounding,
4547             int *flags)
4548{
4549    int expmax = FP64_EXP_BIAS + FP64_BITS - 1;
4550    uint64_t x;
4551    int err;
4552
4553    if (exp > expmax) {
4554        *flags = FPLIB_IOC;
4555        return ((uint64_t)!u << (FP64_BITS - 1)) - !sgn;
4556    }
4557
4558    x = lsr64(mnt << FP64_EXP_BITS, expmax - exp);
4559    err = (exp > expmax - 2 ? 0 :
4560           (lsr64(mnt << FP64_EXP_BITS, expmax - 2 - exp) & 3) |
4561           !!(mnt << FP64_EXP_BITS & (lsl64(1, expmax - 2 - exp) - 1)));
4562
4563    switch (rounding) {
4564      case FPRounding_TIEEVEN:
4565        x += (err == 3 || (err == 2 && (x & 1)));
4566        break;
4567      case FPRounding_POSINF:
4568        x += err && !sgn;
4569        break;
4570      case FPRounding_NEGINF:
4571        x += err && sgn;
4572        break;
4573      case FPRounding_ZERO:
4574        break;
4575      case FPRounding_TIEAWAY:
4576        x += err >> 1;
4577        break;
4578      default:
4579        panic("Unrecognized FP rounding mode");
4580    }
4581
4582    if (u ? sgn && x : x > (1ULL << (FP64_BITS - 1)) - !sgn) {
4583        *flags = FPLIB_IOC;
4584        return ((uint64_t)!u << (FP64_BITS - 1)) - !sgn;
4585    }
4586
4587    if (err) {
4588        *flags = FPLIB_IXC;
4589    }
4590
4591    return sgn ? -x : x;
4592}
4593
4594static uint32_t
4595FPToFixed_32(int sgn, int exp, uint64_t mnt, bool u, FPRounding rounding,
4596             int *flags)
4597{
4598    uint64_t x = FPToFixed_64(sgn, exp, mnt, u, rounding, flags);
4599    if (u ? x >= 1ULL << FP32_BITS :
4600        !(x < 1ULL << (FP32_BITS - 1) ||
4601          (uint64_t)-x <= (uint64_t)1 << (FP32_BITS - 1))) {
4602        *flags = FPLIB_IOC;
4603        x = ((uint32_t)!u << (FP32_BITS - 1)) - !sgn;
4604    }
4605    return x;
4606}
4607
4608static uint16_t
4609FPToFixed_16(int sgn, int exp, uint64_t mnt, bool u, FPRounding rounding,
4610             int *flags)
4611{
4612    uint64_t x = FPToFixed_64(sgn, exp, mnt, u, rounding, flags);
4613    if (u ? x >= 1ULL << FP16_BITS :
4614        !(x < 1ULL << (FP16_BITS - 1) ||
4615          (uint64_t)-x <= (uint64_t)1 << (FP16_BITS - 1))) {
4616        *flags = FPLIB_IOC;
4617        x = ((uint16_t)!u << (FP16_BITS - 1)) - !sgn;
4618    }
4619    return x;
4620}
4621
4622template <>
4623uint16_t
4624fplibFPToFixed(uint16_t op, int fbits, bool u, FPRounding rounding,
4625               FPSCR &fpscr)
4626{
4627    int flags = 0;
4628    int sgn, exp;
4629    uint16_t mnt, result;
4630
4631    // Unpack using FPCR to determine if subnormals are flushed-to-zero:
4632    fp16_unpack(&sgn, &exp, &mnt, op, modeConv(fpscr), &flags);
4633
4634    // If NaN, set cumulative flag or take exception:
4635    if (fp16_is_NaN(exp, mnt)) {
4636        flags = FPLIB_IOC;
4637        result = 0;
4638    } else {
4639        assert(fbits >= 0);
4640        // Infinity is treated as an ordinary normalised number that saturates.
4641        result =
4642            FPToFixed_16(sgn, exp + FP64_EXP_BIAS - FP16_EXP_BIAS + fbits,
4643                         (uint64_t)mnt << (FP64_MANT_BITS - FP16_MANT_BITS),
4644                         u, rounding, &flags);
4645    }
4646
4647    set_fpscr0(fpscr, flags);
4648
4649    return result;
4650}
4651
4652template <>
4653uint32_t
4654fplibFPToFixed(uint16_t op, int fbits, bool u, FPRounding rounding,
4655               FPSCR &fpscr)
4656{
4657    int flags = 0;
4658    int sgn, exp;
4659    uint16_t mnt;
4660    uint32_t result;
4661
4662    // Unpack using FPCR to determine if subnormals are flushed-to-zero:
4663    fp16_unpack(&sgn, &exp, &mnt, op, modeConv(fpscr), &flags);
4664
4665    // If NaN, set cumulative flag or take exception:
4666    if (fp16_is_NaN(exp, mnt)) {
4667        flags = FPLIB_IOC;
4668        result = 0;
4669    } else {
4670        assert(fbits >= 0);
4671        if (exp == FP16_EXP_INF)
4672            exp = 255; // infinity: make it big enough to saturate
4673        result =
4674            FPToFixed_32(sgn, exp + FP64_EXP_BIAS - FP16_EXP_BIAS + fbits,
4675                         (uint64_t)mnt << (FP64_MANT_BITS - FP16_MANT_BITS),
4676                         u, rounding, &flags);
4677    }
4678
4679    set_fpscr0(fpscr, flags);
4680
4681    return result;
4682}
4683
4684template <>
4685uint32_t
4686fplibFPToFixed(uint32_t op, int fbits, bool u, FPRounding rounding, FPSCR &fpscr)
4687{
4688    int flags = 0;
4689    int sgn, exp;
4690    uint32_t mnt, result;
4691
4692    // Unpack using FPCR to determine if subnormals are flushed-to-zero:
4693    fp32_unpack(&sgn, &exp, &mnt, op, modeConv(fpscr), &flags);
4694
4695    // If NaN, set cumulative flag or take exception:
4696    if (fp32_is_NaN(exp, mnt)) {
4697        flags = FPLIB_IOC;
4698        result = 0;
4699    } else {
4700        assert(fbits >= 0);
4701        // Infinity is treated as an ordinary normalised number that saturates.
4702        result =
4703            FPToFixed_32(sgn, exp + FP64_EXP_BIAS - FP32_EXP_BIAS + fbits,
4704                         (uint64_t)mnt << (FP64_MANT_BITS - FP32_MANT_BITS),
4705                         u, rounding, &flags);
4706    }
4707
4708    set_fpscr0(fpscr, flags);
4709
4710    return result;
4711}
4712
4713template <>
4714uint32_t
4715fplibFPToFixed(uint64_t op, int fbits, bool u, FPRounding rounding, FPSCR &fpscr)
4716{
4717    int flags = 0;
4718    int sgn, exp;
4719    uint64_t mnt;
4720    uint32_t result;
4721
4722    // Unpack using FPCR to determine if subnormals are flushed-to-zero:
4723    fp64_unpack(&sgn, &exp, &mnt, op, modeConv(fpscr), &flags);
4724
4725    // If NaN, set cumulative flag or take exception:
4726    if (fp64_is_NaN(exp, mnt)) {
4727        flags = FPLIB_IOC;
4728        result = 0;
4729    } else {
4730        assert(fbits >= 0);
4731        // Infinity is treated as an ordinary normalised number that saturates.
4732        result = FPToFixed_32(sgn, exp + fbits, mnt, u, rounding, &flags);
4733    }
4734
4735    set_fpscr0(fpscr, flags);
4736
4737    return result;
4738}
4739
4740template <>
4741uint64_t
4742fplibFPToFixed(uint16_t op, int fbits, bool u, FPRounding rounding,
4743               FPSCR &fpscr)
4744{
4745    int flags = 0;
4746    int sgn, exp;
4747    uint16_t mnt;
4748    uint64_t result;
4749
4750    // Unpack using FPCR to determine if subnormals are flushed-to-zero:
4751    fp16_unpack(&sgn, &exp, &mnt, op, modeConv(fpscr), &flags);
4752
4753    // If NaN, set cumulative flag or take exception:
4754    if (fp16_is_NaN(exp, mnt)) {
4755        flags = FPLIB_IOC;
4756        result = 0;
4757    } else {
4758        assert(fbits >= 0);
4759        if (exp == FP16_EXP_INF)
4760            exp = 255; // infinity: make it big enough to saturate
4761        result =
4762            FPToFixed_64(sgn, exp + FP64_EXP_BIAS - FP16_EXP_BIAS + fbits,
4763                         (uint64_t)mnt << (FP64_MANT_BITS - FP16_MANT_BITS),
4764                         u, rounding, &flags);
4765    }
4766
4767    set_fpscr0(fpscr, flags);
4768
4769    return result;
4770}
4771
4772template <>
4773uint64_t
4774fplibFPToFixed(uint32_t op, int fbits, bool u, FPRounding rounding, FPSCR &fpscr)
4775{
4776    int flags = 0;
4777    int sgn, exp;
4778    uint32_t mnt;
4779    uint64_t result;
4780
4781    // Unpack using FPCR to determine if subnormals are flushed-to-zero:
4782    fp32_unpack(&sgn, &exp, &mnt, op, modeConv(fpscr), &flags);
4783
4784    // If NaN, set cumulative flag or take exception:
4785    if (fp32_is_NaN(exp, mnt)) {
4786        flags = FPLIB_IOC;
4787        result = 0;
4788    } else {
4789        assert(fbits >= 0);
4790        // Infinity is treated as an ordinary normalised number that saturates.
4791        result =
4792            FPToFixed_64(sgn, exp + FP64_EXP_BIAS - FP32_EXP_BIAS + fbits,
4793                         (uint64_t)mnt << (FP64_MANT_BITS - FP32_MANT_BITS),
4794                         u, rounding, &flags);
4795    }
4796
4797    set_fpscr0(fpscr, flags);
4798
4799    return result;
4800}
4801
4802template <>
4803uint64_t
4804fplibFPToFixed(uint64_t op, int fbits, bool u, FPRounding rounding, FPSCR &fpscr)
4805{
4806    int flags = 0;
4807    int sgn, exp;
4808    uint64_t mnt, result;
4809
4810    // Unpack using FPCR to determine if subnormals are flushed-to-zero:
4811    fp64_unpack(&sgn, &exp, &mnt, op, modeConv(fpscr), &flags);
4812
4813    // If NaN, set cumulative flag or take exception:
4814    if (fp64_is_NaN(exp, mnt)) {
4815        flags = FPLIB_IOC;
4816        result = 0;
4817    } else {
4818        assert(fbits >= 0);
4819        // Infinity is treated as an ordinary normalised number that saturates.
4820        result = FPToFixed_64(sgn, exp + fbits, mnt, u, rounding, &flags);
4821    }
4822
4823    set_fpscr0(fpscr, flags);
4824
4825    return result;
4826}
4827
4828static uint16_t
4829fp16_cvtf(uint64_t a, int fbits, int u, int mode, int *flags)
4830{
4831    int x_sgn = !u && a >> (FP64_BITS - 1);
4832    int x_exp = FP16_EXP_BIAS + FP64_BITS - 1 - fbits;
4833    uint64_t x_mnt = x_sgn ? -a : a;
4834
4835    // Handle zero:
4836    if (!x_mnt) {
4837        return fp16_zero(0);
4838    }
4839
4840    // Normalise into FP16_BITS bits, collapsing error into bottom bit:
4841    x_mnt = fp64_normalise(x_mnt, &x_exp);
4842    x_mnt = (x_mnt >> (FP64_BITS - FP16_BITS - 1) |
4843             !!(x_mnt & ((1ULL << (FP64_BITS - FP16_BITS - 1)) - 1)));
4844
4845    return fp16_round(x_sgn, x_exp, x_mnt, mode, flags);
4846}
4847
4848static uint32_t
4849fp32_cvtf(uint64_t a, int fbits, int u, int mode, int *flags)
4850{
4851    int x_sgn = !u && a >> (FP64_BITS - 1);
4852    int x_exp = FP32_EXP_BIAS + FP64_BITS - 1 - fbits;
4853    uint64_t x_mnt = x_sgn ? -a : a;
4854
4855    // Handle zero:
4856    if (!x_mnt) {
4857        return fp32_zero(0);
4858    }
4859
4860    // Normalise into FP32_BITS bits, collapsing error into bottom bit:
4861    x_mnt = fp64_normalise(x_mnt, &x_exp);
4862    x_mnt = (x_mnt >> (FP64_BITS - FP32_BITS - 1) |
4863             !!(x_mnt & ((1ULL << (FP64_BITS - FP32_BITS - 1)) - 1)));
4864
4865    return fp32_round(x_sgn, x_exp, x_mnt, mode, flags);
4866}
4867
4868static uint64_t
4869fp64_cvtf(uint64_t a, int fbits, int u, int mode, int *flags)
4870{
4871    int x_sgn = !u && a >> (FP64_BITS - 1);
4872    int x_exp = FP64_EXP_BIAS + FP64_BITS - 1 - fbits;
4873    uint64_t x_mnt = x_sgn ? -a : a;
4874
4875    // Handle zero:
4876    if (!x_mnt) {
4877        return fp64_zero(0);
4878    }
4879
4880    x_mnt = fp64_normalise(x_mnt, &x_exp);
4881
4882    return fp64_round(x_sgn, x_exp, x_mnt << 1, mode, flags);
4883}
4884
4885template <>
4886uint16_t
4887fplibFixedToFP(uint64_t op, int fbits, bool u, FPRounding rounding,
4888               FPSCR &fpscr)
4889{
4890    int flags = 0;
4891    uint16_t res = fp16_cvtf(op, fbits, u,
4892                             (int)rounding | ((uint32_t)fpscr >> 22 & 12),
4893                             &flags);
4894    set_fpscr0(fpscr, flags);
4895    return res;
4896}
4897
4898template <>
4899uint32_t
4900fplibFixedToFP(uint64_t op, int fbits, bool u, FPRounding rounding, FPSCR &fpscr)
4901{
4902    int flags = 0;
4903    uint32_t res = fp32_cvtf(op, fbits, u,
4904                             (int)rounding | ((uint32_t)fpscr >> 22 & 12),
4905                             &flags);
4906    set_fpscr0(fpscr, flags);
4907    return res;
4908}
4909
4910template <>
4911uint64_t
4912fplibFixedToFP(uint64_t op, int fbits, bool u, FPRounding rounding, FPSCR &fpscr)
4913{
4914    int flags = 0;
4915    uint64_t res = fp64_cvtf(op, fbits, u,
4916                             (int)rounding | ((uint32_t)fpscr >> 22 & 12),
4917                             &flags);
4918    set_fpscr0(fpscr, flags);
4919    return res;
4920}
4921
4922template <>
4923uint16_t
4924fplibInfinity(int sgn)
4925{
4926    return fp16_infinity(sgn);
4927}
4928
4929template <>
4930uint32_t
4931fplibInfinity(int sgn)
4932{
4933    return fp32_infinity(sgn);
4934}
4935
4936template <>
4937uint64_t
4938fplibInfinity(int sgn)
4939{
4940    return fp64_infinity(sgn);
4941}
4942
4943template <>
4944uint16_t
4945fplibDefaultNaN()
4946{
4947    return fp16_defaultNaN();
4948}
4949
4950template <>
4951uint32_t
4952fplibDefaultNaN()
4953{
4954    return fp32_defaultNaN();
4955}
4956
4957template <>
4958uint64_t
4959fplibDefaultNaN()
4960{
4961    return fp64_defaultNaN();
4962}
4963
4964}
4965