3 fp_arith.c: floating-point math routines for the Linux-m68k
4 floating point emulator.
6 Copyright (c) 1998-1999 David Huggins-Daines.
8 Somewhat based on the AlphaLinux floating point emulator, by David
11 You may copy, modify, and redistribute this file under the terms of
12 the GNU General Public License, version 2, or any later version, at
17 #include "multi_arith.h"
20 const struct fp_ext fp_QNaN =
25 const struct fp_ext fp_Inf =
30 /* let's start with the easy ones */
33 fp_fabs(struct fp_ext *dest, struct fp_ext *src)
35 dprint(PINSTR, "fabs\n");
37 fp_monadic_check(dest, src);
45 fp_fneg(struct fp_ext *dest, struct fp_ext *src)
47 dprint(PINSTR, "fneg\n");
49 fp_monadic_check(dest, src);
51 dest->sign = !dest->sign;
56 /* Now, the slightly harder ones */
58 /* fp_fadd: Implements the kernel of the FADD, FSADD, FDADD, FSUB,
59 FDSUB, and FCMP instructions. */
62 fp_fadd(struct fp_ext *dest, struct fp_ext *src)
66 dprint(PINSTR, "fadd\n");
68 fp_dyadic_check(dest, src);
71 /* infinity - infinity == NaN */
72 if (IS_INF(src) && (src->sign != dest->sign))
77 fp_copy_ext(dest, src);
83 if (src->sign != dest->sign) {
84 if (FPDATA->rnd == FPCR_ROUND_RM)
90 fp_copy_ext(dest, src);
94 dest->lowmant = src->lowmant = 0;
96 if ((diff = dest->exp - src->exp) > 0)
97 fp_denormalize(src, diff);
98 else if ((diff = -diff) > 0)
99 fp_denormalize(dest, diff);
101 if (dest->sign == src->sign) {
102 if (fp_addmant(dest, src))
103 if (!fp_addcarry(dest))
106 if (dest->mant.m64 < src->mant.m64) {
107 fp_submant(dest, src, dest);
108 dest->sign = !dest->sign;
110 fp_submant(dest, dest, src);
116 /* fp_fsub: Implements the kernel of the FSUB, FSSUB, and FDSUB
119 Remember that the arguments are in assembler-syntax order! */
122 fp_fsub(struct fp_ext *dest, struct fp_ext *src)
124 dprint(PINSTR, "fsub ");
126 src->sign = !src->sign;
127 return fp_fadd(dest, src);
132 fp_fcmp(struct fp_ext *dest, struct fp_ext *src)
134 dprint(PINSTR, "fcmp ");
136 FPDATA->temp[1] = *dest;
137 src->sign = !src->sign;
138 return fp_fadd(&FPDATA->temp[1], src);
142 fp_ftst(struct fp_ext *dest, struct fp_ext *src)
144 dprint(PINSTR, "ftst\n");
152 fp_fmul(struct fp_ext *dest, struct fp_ext *src)
154 union fp_mant128 temp;
157 dprint(PINSTR, "fmul\n");
159 fp_dyadic_check(dest, src);
161 /* calculate the correct sign now, as it's necessary for infinities */
162 dest->sign = src->sign ^ dest->sign;
164 /* Handle infinities */
174 fp_copy_ext(dest, src);
178 /* Of course, as we all know, zero * anything = zero. You may
179 not have known that it might be a positive or negative
181 if (IS_ZERO(dest) || IS_ZERO(src)) {
189 exp = dest->exp + src->exp - 0x3ffe;
191 /* shift up the mantissa for denormalized numbers,
192 so that the highest bit is set, this makes the
193 shift of the result below easier */
194 if ((long)dest->mant.m32[0] >= 0)
195 exp -= fp_overnormalize(dest);
196 if ((long)src->mant.m32[0] >= 0)
197 exp -= fp_overnormalize(src);
199 /* now, do a 64-bit multiply with expansion */
200 fp_multiplymant(&temp, dest, src);
202 /* normalize it back to 64 bits and stuff it back into the
203 destination struct */
204 if ((long)temp.m32[0] > 0) {
206 fp_putmant128(dest, &temp, 1);
208 fp_putmant128(dest, &temp, 0);
216 fp_set_sr(FPSR_EXC_UNFL);
217 fp_denormalize(dest, -exp);
223 /* fp_fdiv: Implements the "kernel" of the FDIV, FSDIV, FDDIV and
224 FSGLDIV instructions.
226 Note that the order of the operands is counter-intuitive: instead
227 of src / dest, the result is actually dest / src. */
230 fp_fdiv(struct fp_ext *dest, struct fp_ext *src)
232 union fp_mant128 temp;
235 dprint(PINSTR, "fdiv\n");
237 fp_dyadic_check(dest, src);
239 /* calculate the correct sign now, as it's necessary for infinities */
240 dest->sign = src->sign ^ dest->sign;
242 /* Handle infinities */
244 /* infinity / infinity = NaN (quiet, as always) */
247 /* infinity / anything else = infinity (with approprate sign) */
251 /* anything / infinity = zero (with appropriate sign) */
261 /* zero / zero = NaN */
264 /* zero / anything else = zero */
268 /* anything / zero = infinity (with appropriate sign) */
269 fp_set_sr(FPSR_EXC_DZ);
276 exp = dest->exp - src->exp + 0x3fff;
278 /* shift up the mantissa for denormalized numbers,
279 so that the highest bit is set, this makes lots
280 of things below easier */
281 if ((long)dest->mant.m32[0] >= 0)
282 exp -= fp_overnormalize(dest);
283 if ((long)src->mant.m32[0] >= 0)
284 exp -= fp_overnormalize(src);
286 /* now, do the 64-bit divide */
287 fp_dividemant(&temp, dest, src);
289 /* normalize it back to 64 bits and stuff it back into the
290 destination struct */
293 fp_putmant128(dest, &temp, 32);
295 fp_putmant128(dest, &temp, 31);
303 fp_set_sr(FPSR_EXC_UNFL);
304 fp_denormalize(dest, -exp);
311 fp_fsglmul(struct fp_ext *dest, struct fp_ext *src)
315 dprint(PINSTR, "fsglmul\n");
317 fp_dyadic_check(dest, src);
319 /* calculate the correct sign now, as it's necessary for infinities */
320 dest->sign = src->sign ^ dest->sign;
322 /* Handle infinities */
332 fp_copy_ext(dest, src);
336 /* Of course, as we all know, zero * anything = zero. You may
337 not have known that it might be a positive or negative
339 if (IS_ZERO(dest) || IS_ZERO(src)) {
347 exp = dest->exp + src->exp - 0x3ffe;
349 /* do a 32-bit multiply */
350 fp_mul64(dest->mant.m32[0], dest->mant.m32[1],
351 dest->mant.m32[0] & 0xffffff00,
352 src->mant.m32[0] & 0xffffff00);
360 fp_set_sr(FPSR_EXC_UNFL);
361 fp_denormalize(dest, -exp);
368 fp_fsgldiv(struct fp_ext *dest, struct fp_ext *src)
371 unsigned long quot, rem;
373 dprint(PINSTR, "fsgldiv\n");
375 fp_dyadic_check(dest, src);
377 /* calculate the correct sign now, as it's necessary for infinities */
378 dest->sign = src->sign ^ dest->sign;
380 /* Handle infinities */
382 /* infinity / infinity = NaN (quiet, as always) */
385 /* infinity / anything else = infinity (with approprate sign) */
389 /* anything / infinity = zero (with appropriate sign) */
399 /* zero / zero = NaN */
402 /* zero / anything else = zero */
406 /* anything / zero = infinity (with appropriate sign) */
407 fp_set_sr(FPSR_EXC_DZ);
414 exp = dest->exp - src->exp + 0x3fff;
416 dest->mant.m32[0] &= 0xffffff00;
417 src->mant.m32[0] &= 0xffffff00;
419 /* do the 32-bit divide */
420 if (dest->mant.m32[0] >= src->mant.m32[0]) {
421 fp_sub64(dest->mant, src->mant);
422 fp_div64(quot, rem, dest->mant.m32[0], 0, src->mant.m32[0]);
423 dest->mant.m32[0] = 0x80000000 | (quot >> 1);
424 dest->mant.m32[1] = (quot & 1) | rem; /* only for rounding */
426 fp_div64(quot, rem, dest->mant.m32[0], 0, src->mant.m32[0]);
427 dest->mant.m32[0] = quot;
428 dest->mant.m32[1] = rem; /* only for rounding */
438 fp_set_sr(FPSR_EXC_UNFL);
439 fp_denormalize(dest, -exp);
445 /* fp_roundint: Internal rounding function for use by several of these
446 emulated instructions.
448 This one rounds off the fractional part using the rounding mode
451 static void fp_roundint(struct fp_ext *dest, int mode)
453 union fp_mant64 oldmant;
456 if (!fp_normalize_ext(dest))
459 /* infinities and zeroes */
460 if (IS_INF(dest) || IS_ZERO(dest))
463 /* first truncate the lower bits */
464 oldmant = dest->mant;
469 case 0x3fff ... 0x401e:
470 dest->mant.m32[0] &= 0xffffffffU << (0x401e - dest->exp);
471 dest->mant.m32[1] = 0;
472 if (oldmant.m64 == dest->mant.m64)
475 case 0x401f ... 0x403e:
476 dest->mant.m32[1] &= 0xffffffffU << (0x403e - dest->exp);
477 if (oldmant.m32[1] == dest->mant.m32[1])
483 fp_set_sr(FPSR_EXC_INEX2);
485 /* We might want to normalize upwards here... however, since
486 we know that this is only called on the output of fp_fdiv,
487 or with the input to fp_fint or fp_fintrz, and the inputs
488 to all these functions are either normal or denormalized
489 (no subnormals allowed!), there's really no need.
491 In the case of fp_fdiv, observe that 0x80000000 / 0xffff =
492 0xffff8000, and the same holds for 128-bit / 64-bit. (i.e. the
493 smallest possible normal dividend and the largest possible normal
494 divisor will still produce a normal quotient, therefore, (normal
495 << 64) / normal is normal in all cases) */
503 /* As noted above, the input is always normal, so the
504 guard bit (bit 63) is always set. therefore, the
505 only case in which we will NOT round to 1.0 is when
506 the input is exactly 0.5. */
507 if (oldmant.m64 == (1ULL << 63))
510 case 0x3fff ... 0x401d:
511 mask = 1 << (0x401d - dest->exp);
512 if (!(oldmant.m32[0] & mask))
514 if (oldmant.m32[0] & (mask << 1))
516 if (!(oldmant.m32[0] << (dest->exp - 0x3ffd)) &&
521 if (!(oldmant.m32[1] >= 0))
523 if (oldmant.m32[0] & 1)
525 if (!(oldmant.m32[1] << 1))
528 case 0x401f ... 0x403d:
529 mask = 1 << (0x403d - dest->exp);
530 if (!(oldmant.m32[1] & mask))
532 if (oldmant.m32[1] & (mask << 1))
534 if (!(oldmant.m32[1] << (dest->exp - 0x401d)))
544 if (dest->sign ^ (mode - FPCR_ROUND_RM))
552 dest->mant.m64 = 1ULL << 63;
554 case 0x3fff ... 0x401e:
555 mask = 1 << (0x401e - dest->exp);
556 if (dest->mant.m32[0] += mask)
558 dest->mant.m32[0] = 0x80000000;
561 case 0x401f ... 0x403e:
562 mask = 1 << (0x403e - dest->exp);
563 if (dest->mant.m32[1] += mask)
565 if (dest->mant.m32[0] += 1)
567 dest->mant.m32[0] = 0x80000000;
573 /* modrem_kernel: Implementation of the FREM and FMOD instructions
574 (which are exactly the same, except for the rounding used on the
575 intermediate value) */
577 static struct fp_ext *
578 modrem_kernel(struct fp_ext *dest, struct fp_ext *src, int mode)
582 fp_dyadic_check(dest, src);
584 /* Infinities and zeros */
585 if (IS_INF(dest) || IS_ZERO(src)) {
589 if (IS_ZERO(dest) || IS_INF(src))
592 /* FIXME: there is almost certainly a smarter way to do this */
593 fp_copy_ext(&tmp, dest);
594 fp_fdiv(&tmp, src); /* NOTE: src might be modified */
595 fp_roundint(&tmp, mode);
599 /* set the quotient byte */
600 fp_set_quotient((dest->mant.m64 & 0x7f) | (dest->sign << 7));
604 /* fp_fmod: Implements the kernel of the FMOD instruction.
606 Again, the argument order is backwards. The result, as defined in
607 the Motorola manuals, is:
609 fmod(src,dest) = (dest - (src * floor(dest / src))) */
612 fp_fmod(struct fp_ext *dest, struct fp_ext *src)
614 dprint(PINSTR, "fmod\n");
615 return modrem_kernel(dest, src, FPCR_ROUND_RZ);
618 /* fp_frem: Implements the kernel of the FREM instruction.
620 frem(src,dest) = (dest - (src * round(dest / src)))
624 fp_frem(struct fp_ext *dest, struct fp_ext *src)
626 dprint(PINSTR, "frem\n");
627 return modrem_kernel(dest, src, FPCR_ROUND_RN);
631 fp_fint(struct fp_ext *dest, struct fp_ext *src)
633 dprint(PINSTR, "fint\n");
635 fp_copy_ext(dest, src);
637 fp_roundint(dest, FPDATA->rnd);
643 fp_fintrz(struct fp_ext *dest, struct fp_ext *src)
645 dprint(PINSTR, "fintrz\n");
647 fp_copy_ext(dest, src);
649 fp_roundint(dest, FPCR_ROUND_RZ);
655 fp_fscale(struct fp_ext *dest, struct fp_ext *src)
659 dprint(PINSTR, "fscale\n");
661 fp_dyadic_check(dest, src);
672 if (IS_ZERO(src) || IS_ZERO(dest))
675 /* Source exponent out of range */
676 if (src->exp >= 0x400c) {
681 /* src must be rounded with round to zero. */
682 oldround = FPDATA->rnd;
683 FPDATA->rnd = FPCR_ROUND_RZ;
684 scale = fp_conv_ext2long(src);
685 FPDATA->rnd = oldround;
690 if (scale >= 0x7fff) {
692 } else if (scale <= 0) {
693 fp_set_sr(FPSR_EXC_UNFL);
694 fp_denormalize(dest, -scale);