1 // SPDX-License-Identifier: GPL-2.0-or-later
4 fp_arith.c: floating-point math routines for the Linux-m68k
5 floating point emulator.
7 Copyright (c) 1998-1999 David Huggins-Daines.
9 Somewhat based on the AlphaLinux floating point emulator, by David
15 #include "multi_arith.h"
18 const struct fp_ext fp_QNaN
=
24 const struct fp_ext fp_Inf
=
29 /* let's start with the easy ones */
32 fp_fabs(struct fp_ext
*dest
, struct fp_ext
*src
)
34 dprint(PINSTR
, "fabs\n");
36 fp_monadic_check(dest
, src
);
44 fp_fneg(struct fp_ext
*dest
, struct fp_ext
*src
)
46 dprint(PINSTR
, "fneg\n");
48 fp_monadic_check(dest
, src
);
50 dest
->sign
= !dest
->sign
;
55 /* Now, the slightly harder ones */
57 /* fp_fadd: Implements the kernel of the FADD, FSADD, FDADD, FSUB,
58 FDSUB, and FCMP instructions. */
61 fp_fadd(struct fp_ext
*dest
, struct fp_ext
*src
)
65 dprint(PINSTR
, "fadd\n");
67 fp_dyadic_check(dest
, src
);
70 /* infinity - infinity == NaN */
71 if (IS_INF(src
) && (src
->sign
!= dest
->sign
))
76 fp_copy_ext(dest
, src
);
82 if (src
->sign
!= dest
->sign
) {
83 if (FPDATA
->rnd
== FPCR_ROUND_RM
)
89 fp_copy_ext(dest
, src
);
93 dest
->lowmant
= src
->lowmant
= 0;
95 if ((diff
= dest
->exp
- src
->exp
) > 0)
96 fp_denormalize(src
, diff
);
97 else if ((diff
= -diff
) > 0)
98 fp_denormalize(dest
, diff
);
100 if (dest
->sign
== src
->sign
) {
101 if (fp_addmant(dest
, src
))
102 if (!fp_addcarry(dest
))
105 if (dest
->mant
.m64
< src
->mant
.m64
) {
106 fp_submant(dest
, src
, dest
);
107 dest
->sign
= !dest
->sign
;
109 fp_submant(dest
, dest
, src
);
115 /* fp_fsub: Implements the kernel of the FSUB, FSSUB, and FDSUB
118 Remember that the arguments are in assembler-syntax order! */
121 fp_fsub(struct fp_ext
*dest
, struct fp_ext
*src
)
123 dprint(PINSTR
, "fsub ");
125 src
->sign
= !src
->sign
;
126 return fp_fadd(dest
, src
);
131 fp_fcmp(struct fp_ext
*dest
, struct fp_ext
*src
)
133 dprint(PINSTR
, "fcmp ");
135 FPDATA
->temp
[1] = *dest
;
136 src
->sign
= !src
->sign
;
137 return fp_fadd(&FPDATA
->temp
[1], src
);
141 fp_ftst(struct fp_ext
*dest
, struct fp_ext
*src
)
143 dprint(PINSTR
, "ftst\n");
151 fp_fmul(struct fp_ext
*dest
, struct fp_ext
*src
)
153 union fp_mant128 temp
;
156 dprint(PINSTR
, "fmul\n");
158 fp_dyadic_check(dest
, src
);
160 /* calculate the correct sign now, as it's necessary for infinities */
161 dest
->sign
= src
->sign
^ dest
->sign
;
163 /* Handle infinities */
173 fp_copy_ext(dest
, src
);
177 /* Of course, as we all know, zero * anything = zero. You may
178 not have known that it might be a positive or negative
180 if (IS_ZERO(dest
) || IS_ZERO(src
)) {
188 exp
= dest
->exp
+ src
->exp
- 0x3ffe;
190 /* shift up the mantissa for denormalized numbers,
191 so that the highest bit is set, this makes the
192 shift of the result below easier */
193 if ((long)dest
->mant
.m32
[0] >= 0)
194 exp
-= fp_overnormalize(dest
);
195 if ((long)src
->mant
.m32
[0] >= 0)
196 exp
-= fp_overnormalize(src
);
198 /* now, do a 64-bit multiply with expansion */
199 fp_multiplymant(&temp
, dest
, src
);
201 /* normalize it back to 64 bits and stuff it back into the
202 destination struct */
203 if ((long)temp
.m32
[0] > 0) {
205 fp_putmant128(dest
, &temp
, 1);
207 fp_putmant128(dest
, &temp
, 0);
215 fp_set_sr(FPSR_EXC_UNFL
);
216 fp_denormalize(dest
, -exp
);
222 /* fp_fdiv: Implements the "kernel" of the FDIV, FSDIV, FDDIV and
223 FSGLDIV instructions.
225 Note that the order of the operands is counter-intuitive: instead
226 of src / dest, the result is actually dest / src. */
229 fp_fdiv(struct fp_ext
*dest
, struct fp_ext
*src
)
231 union fp_mant128 temp
;
234 dprint(PINSTR
, "fdiv\n");
236 fp_dyadic_check(dest
, src
);
238 /* calculate the correct sign now, as it's necessary for infinities */
239 dest
->sign
= src
->sign
^ dest
->sign
;
241 /* Handle infinities */
243 /* infinity / infinity = NaN (quiet, as always) */
246 /* infinity / anything else = infinity (with approprate sign) */
250 /* anything / infinity = zero (with appropriate sign) */
260 /* zero / zero = NaN */
263 /* zero / anything else = zero */
267 /* anything / zero = infinity (with appropriate sign) */
268 fp_set_sr(FPSR_EXC_DZ
);
275 exp
= dest
->exp
- src
->exp
+ 0x3fff;
277 /* shift up the mantissa for denormalized numbers,
278 so that the highest bit is set, this makes lots
279 of things below easier */
280 if ((long)dest
->mant
.m32
[0] >= 0)
281 exp
-= fp_overnormalize(dest
);
282 if ((long)src
->mant
.m32
[0] >= 0)
283 exp
-= fp_overnormalize(src
);
285 /* now, do the 64-bit divide */
286 fp_dividemant(&temp
, dest
, src
);
288 /* normalize it back to 64 bits and stuff it back into the
289 destination struct */
292 fp_putmant128(dest
, &temp
, 32);
294 fp_putmant128(dest
, &temp
, 31);
302 fp_set_sr(FPSR_EXC_UNFL
);
303 fp_denormalize(dest
, -exp
);
310 fp_fsglmul(struct fp_ext
*dest
, struct fp_ext
*src
)
314 dprint(PINSTR
, "fsglmul\n");
316 fp_dyadic_check(dest
, src
);
318 /* calculate the correct sign now, as it's necessary for infinities */
319 dest
->sign
= src
->sign
^ dest
->sign
;
321 /* Handle infinities */
331 fp_copy_ext(dest
, src
);
335 /* Of course, as we all know, zero * anything = zero. You may
336 not have known that it might be a positive or negative
338 if (IS_ZERO(dest
) || IS_ZERO(src
)) {
346 exp
= dest
->exp
+ src
->exp
- 0x3ffe;
348 /* do a 32-bit multiply */
349 fp_mul64(dest
->mant
.m32
[0], dest
->mant
.m32
[1],
350 dest
->mant
.m32
[0] & 0xffffff00,
351 src
->mant
.m32
[0] & 0xffffff00);
359 fp_set_sr(FPSR_EXC_UNFL
);
360 fp_denormalize(dest
, -exp
);
367 fp_fsgldiv(struct fp_ext
*dest
, struct fp_ext
*src
)
370 unsigned long quot
, rem
;
372 dprint(PINSTR
, "fsgldiv\n");
374 fp_dyadic_check(dest
, src
);
376 /* calculate the correct sign now, as it's necessary for infinities */
377 dest
->sign
= src
->sign
^ dest
->sign
;
379 /* Handle infinities */
381 /* infinity / infinity = NaN (quiet, as always) */
384 /* infinity / anything else = infinity (with approprate sign) */
388 /* anything / infinity = zero (with appropriate sign) */
398 /* zero / zero = NaN */
401 /* zero / anything else = zero */
405 /* anything / zero = infinity (with appropriate sign) */
406 fp_set_sr(FPSR_EXC_DZ
);
413 exp
= dest
->exp
- src
->exp
+ 0x3fff;
415 dest
->mant
.m32
[0] &= 0xffffff00;
416 src
->mant
.m32
[0] &= 0xffffff00;
418 /* do the 32-bit divide */
419 if (dest
->mant
.m32
[0] >= src
->mant
.m32
[0]) {
420 fp_sub64(dest
->mant
, src
->mant
);
421 fp_div64(quot
, rem
, dest
->mant
.m32
[0], 0, src
->mant
.m32
[0]);
422 dest
->mant
.m32
[0] = 0x80000000 | (quot
>> 1);
423 dest
->mant
.m32
[1] = (quot
& 1) | rem
; /* only for rounding */
425 fp_div64(quot
, rem
, dest
->mant
.m32
[0], 0, src
->mant
.m32
[0]);
426 dest
->mant
.m32
[0] = quot
;
427 dest
->mant
.m32
[1] = rem
; /* only for rounding */
437 fp_set_sr(FPSR_EXC_UNFL
);
438 fp_denormalize(dest
, -exp
);
444 /* fp_roundint: Internal rounding function for use by several of these
445 emulated instructions.
447 This one rounds off the fractional part using the rounding mode
450 static void fp_roundint(struct fp_ext
*dest
, int mode
)
452 union fp_mant64 oldmant
;
455 if (!fp_normalize_ext(dest
))
458 /* infinities and zeroes */
459 if (IS_INF(dest
) || IS_ZERO(dest
))
462 /* first truncate the lower bits */
463 oldmant
= dest
->mant
;
468 case 0x3fff ... 0x401e:
469 dest
->mant
.m32
[0] &= 0xffffffffU
<< (0x401e - dest
->exp
);
470 dest
->mant
.m32
[1] = 0;
471 if (oldmant
.m64
== dest
->mant
.m64
)
474 case 0x401f ... 0x403e:
475 dest
->mant
.m32
[1] &= 0xffffffffU
<< (0x403e - dest
->exp
);
476 if (oldmant
.m32
[1] == dest
->mant
.m32
[1])
482 fp_set_sr(FPSR_EXC_INEX2
);
484 /* We might want to normalize upwards here... however, since
485 we know that this is only called on the output of fp_fdiv,
486 or with the input to fp_fint or fp_fintrz, and the inputs
487 to all these functions are either normal or denormalized
488 (no subnormals allowed!), there's really no need.
490 In the case of fp_fdiv, observe that 0x80000000 / 0xffff =
491 0xffff8000, and the same holds for 128-bit / 64-bit. (i.e. the
492 smallest possible normal dividend and the largest possible normal
493 divisor will still produce a normal quotient, therefore, (normal
494 << 64) / normal is normal in all cases) */
502 /* As noted above, the input is always normal, so the
503 guard bit (bit 63) is always set. therefore, the
504 only case in which we will NOT round to 1.0 is when
505 the input is exactly 0.5. */
506 if (oldmant
.m64
== (1ULL << 63))
509 case 0x3fff ... 0x401d:
510 mask
= 1 << (0x401d - dest
->exp
);
511 if (!(oldmant
.m32
[0] & mask
))
513 if (oldmant
.m32
[0] & (mask
<< 1))
515 if (!(oldmant
.m32
[0] << (dest
->exp
- 0x3ffd)) &&
520 if (oldmant
.m32
[1] & 0x80000000)
522 if (oldmant
.m32
[0] & 1)
524 if (!(oldmant
.m32
[1] << 1))
527 case 0x401f ... 0x403d:
528 mask
= 1 << (0x403d - dest
->exp
);
529 if (!(oldmant
.m32
[1] & mask
))
531 if (oldmant
.m32
[1] & (mask
<< 1))
533 if (!(oldmant
.m32
[1] << (dest
->exp
- 0x401d)))
543 if (dest
->sign
^ (mode
- FPCR_ROUND_RM
))
551 dest
->mant
.m64
= 1ULL << 63;
553 case 0x3fff ... 0x401e:
554 mask
= 1 << (0x401e - dest
->exp
);
555 if (dest
->mant
.m32
[0] += mask
)
557 dest
->mant
.m32
[0] = 0x80000000;
560 case 0x401f ... 0x403e:
561 mask
= 1 << (0x403e - dest
->exp
);
562 if (dest
->mant
.m32
[1] += mask
)
564 if (dest
->mant
.m32
[0] += 1)
566 dest
->mant
.m32
[0] = 0x80000000;
572 /* modrem_kernel: Implementation of the FREM and FMOD instructions
573 (which are exactly the same, except for the rounding used on the
574 intermediate value) */
576 static struct fp_ext
*
577 modrem_kernel(struct fp_ext
*dest
, struct fp_ext
*src
, int mode
)
581 fp_dyadic_check(dest
, src
);
583 /* Infinities and zeros */
584 if (IS_INF(dest
) || IS_ZERO(src
)) {
588 if (IS_ZERO(dest
) || IS_INF(src
))
591 /* FIXME: there is almost certainly a smarter way to do this */
592 fp_copy_ext(&tmp
, dest
);
593 fp_fdiv(&tmp
, src
); /* NOTE: src might be modified */
594 fp_roundint(&tmp
, mode
);
598 /* set the quotient byte */
599 fp_set_quotient((dest
->mant
.m64
& 0x7f) | (dest
->sign
<< 7));
603 /* fp_fmod: Implements the kernel of the FMOD instruction.
605 Again, the argument order is backwards. The result, as defined in
606 the Motorola manuals, is:
608 fmod(src,dest) = (dest - (src * floor(dest / src))) */
611 fp_fmod(struct fp_ext
*dest
, struct fp_ext
*src
)
613 dprint(PINSTR
, "fmod\n");
614 return modrem_kernel(dest
, src
, FPCR_ROUND_RZ
);
617 /* fp_frem: Implements the kernel of the FREM instruction.
619 frem(src,dest) = (dest - (src * round(dest / src)))
623 fp_frem(struct fp_ext
*dest
, struct fp_ext
*src
)
625 dprint(PINSTR
, "frem\n");
626 return modrem_kernel(dest
, src
, FPCR_ROUND_RN
);
630 fp_fint(struct fp_ext
*dest
, struct fp_ext
*src
)
632 dprint(PINSTR
, "fint\n");
634 fp_copy_ext(dest
, src
);
636 fp_roundint(dest
, FPDATA
->rnd
);
642 fp_fintrz(struct fp_ext
*dest
, struct fp_ext
*src
)
644 dprint(PINSTR
, "fintrz\n");
646 fp_copy_ext(dest
, src
);
648 fp_roundint(dest
, FPCR_ROUND_RZ
);
654 fp_fscale(struct fp_ext
*dest
, struct fp_ext
*src
)
658 dprint(PINSTR
, "fscale\n");
660 fp_dyadic_check(dest
, src
);
671 if (IS_ZERO(src
) || IS_ZERO(dest
))
674 /* Source exponent out of range */
675 if (src
->exp
>= 0x400c) {
680 /* src must be rounded with round to zero. */
681 oldround
= FPDATA
->rnd
;
682 FPDATA
->rnd
= FPCR_ROUND_RZ
;
683 scale
= fp_conv_ext2long(src
);
684 FPDATA
->rnd
= oldround
;
689 if (scale
>= 0x7fff) {
691 } else if (scale
<= 0) {
692 fp_set_sr(FPSR_EXC_UNFL
);
693 fp_denormalize(dest
, -scale
);