1 /*
2
3 fp_arith.c: floating-point math routines for the Linux-m68k
4 floating point emulator.
5
6 Copyright (c) 1998-1999 David Huggins-Daines.
7
8 Somewhat based on the AlphaLinux floating point emulator, by David
9 Mosberger-Tang.
10
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
13 your convenience.
14 */
15
16 #include "fp_emu.h"
17 #include "multi_arith.h"
18 #include "fp_arith.h"
19
20 const struct fp_ext fp_QNaN =
21 {
22 0, 0, 0x7fff, { ~0 }
23 };
24
25 const struct fp_ext fp_Inf =
26 {
27 0, 0, 0x7fff, { 0 }
28 };
29
30 /* let's start with the easy ones */
31
32 struct fp_ext *
fp_fabs(struct fp_ext * dest,struct fp_ext * src)33 fp_fabs(struct fp_ext *dest, struct fp_ext *src)
34 {
35 dprint(PINSTR, "fabs\n");
36
37 fp_monadic_check(dest, src);
38
39 dest->sign = 0;
40
41 return dest;
42 }
43
44 struct fp_ext *
fp_fneg(struct fp_ext * dest,struct fp_ext * src)45 fp_fneg(struct fp_ext *dest, struct fp_ext *src)
46 {
47 dprint(PINSTR, "fneg\n");
48
49 fp_monadic_check(dest, src);
50
51 dest->sign = !dest->sign;
52
53 return dest;
54 }
55
56 /* Now, the slightly harder ones */
57
58 /* fp_fadd: Implements the kernel of the FADD, FSADD, FDADD, FSUB,
59 FDSUB, and FCMP instructions. */
60
61 struct fp_ext *
fp_fadd(struct fp_ext * dest,struct fp_ext * src)62 fp_fadd(struct fp_ext *dest, struct fp_ext *src)
63 {
64 int diff;
65
66 dprint(PINSTR, "fadd\n");
67
68 fp_dyadic_check(dest, src);
69
70 if (IS_INF(dest)) {
71 /* infinity - infinity == NaN */
72 if (IS_INF(src) && (src->sign != dest->sign))
73 fp_set_nan(dest);
74 return dest;
75 }
76 if (IS_INF(src)) {
77 fp_copy_ext(dest, src);
78 return dest;
79 }
80
81 if (IS_ZERO(dest)) {
82 if (IS_ZERO(src)) {
83 if (src->sign != dest->sign) {
84 if (FPDATA->rnd == FPCR_ROUND_RM)
85 dest->sign = 1;
86 else
87 dest->sign = 0;
88 }
89 } else
90 fp_copy_ext(dest, src);
91 return dest;
92 }
93
94 dest->lowmant = src->lowmant = 0;
95
96 if ((diff = dest->exp - src->exp) > 0)
97 fp_denormalize(src, diff);
98 else if ((diff = -diff) > 0)
99 fp_denormalize(dest, diff);
100
101 if (dest->sign == src->sign) {
102 if (fp_addmant(dest, src))
103 if (!fp_addcarry(dest))
104 return dest;
105 } else {
106 if (dest->mant.m64 < src->mant.m64) {
107 fp_submant(dest, src, dest);
108 dest->sign = !dest->sign;
109 } else
110 fp_submant(dest, dest, src);
111 }
112
113 return dest;
114 }
115
116 /* fp_fsub: Implements the kernel of the FSUB, FSSUB, and FDSUB
117 instructions.
118
119 Remember that the arguments are in assembler-syntax order! */
120
121 struct fp_ext *
fp_fsub(struct fp_ext * dest,struct fp_ext * src)122 fp_fsub(struct fp_ext *dest, struct fp_ext *src)
123 {
124 dprint(PINSTR, "fsub ");
125
126 src->sign = !src->sign;
127 return fp_fadd(dest, src);
128 }
129
130
131 struct fp_ext *
fp_fcmp(struct fp_ext * dest,struct fp_ext * src)132 fp_fcmp(struct fp_ext *dest, struct fp_ext *src)
133 {
134 dprint(PINSTR, "fcmp ");
135
136 FPDATA->temp[1] = *dest;
137 src->sign = !src->sign;
138 return fp_fadd(&FPDATA->temp[1], src);
139 }
140
141 struct fp_ext *
fp_ftst(struct fp_ext * dest,struct fp_ext * src)142 fp_ftst(struct fp_ext *dest, struct fp_ext *src)
143 {
144 dprint(PINSTR, "ftst\n");
145
146 (void)dest;
147
148 return src;
149 }
150
151 struct fp_ext *
fp_fmul(struct fp_ext * dest,struct fp_ext * src)152 fp_fmul(struct fp_ext *dest, struct fp_ext *src)
153 {
154 union fp_mant128 temp;
155 int exp;
156
157 dprint(PINSTR, "fmul\n");
158
159 fp_dyadic_check(dest, src);
160
161 /* calculate the correct sign now, as it's necessary for infinities */
162 dest->sign = src->sign ^ dest->sign;
163
164 /* Handle infinities */
165 if (IS_INF(dest)) {
166 if (IS_ZERO(src))
167 fp_set_nan(dest);
168 return dest;
169 }
170 if (IS_INF(src)) {
171 if (IS_ZERO(dest))
172 fp_set_nan(dest);
173 else
174 fp_copy_ext(dest, src);
175 return dest;
176 }
177
178 /* Of course, as we all know, zero * anything = zero. You may
179 not have known that it might be a positive or negative
180 zero... */
181 if (IS_ZERO(dest) || IS_ZERO(src)) {
182 dest->exp = 0;
183 dest->mant.m64 = 0;
184 dest->lowmant = 0;
185
186 return dest;
187 }
188
189 exp = dest->exp + src->exp - 0x3ffe;
190
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);
198
199 /* now, do a 64-bit multiply with expansion */
200 fp_multiplymant(&temp, dest, src);
201
202 /* normalize it back to 64 bits and stuff it back into the
203 destination struct */
204 if ((long)temp.m32[0] > 0) {
205 exp--;
206 fp_putmant128(dest, &temp, 1);
207 } else
208 fp_putmant128(dest, &temp, 0);
209
210 if (exp >= 0x7fff) {
211 fp_set_ovrflw(dest);
212 return dest;
213 }
214 dest->exp = exp;
215 if (exp < 0) {
216 fp_set_sr(FPSR_EXC_UNFL);
217 fp_denormalize(dest, -exp);
218 }
219
220 return dest;
221 }
222
223 /* fp_fdiv: Implements the "kernel" of the FDIV, FSDIV, FDDIV and
224 FSGLDIV instructions.
225
226 Note that the order of the operands is counter-intuitive: instead
227 of src / dest, the result is actually dest / src. */
228
229 struct fp_ext *
fp_fdiv(struct fp_ext * dest,struct fp_ext * src)230 fp_fdiv(struct fp_ext *dest, struct fp_ext *src)
231 {
232 union fp_mant128 temp;
233 int exp;
234
235 dprint(PINSTR, "fdiv\n");
236
237 fp_dyadic_check(dest, src);
238
239 /* calculate the correct sign now, as it's necessary for infinities */
240 dest->sign = src->sign ^ dest->sign;
241
242 /* Handle infinities */
243 if (IS_INF(dest)) {
244 /* infinity / infinity = NaN (quiet, as always) */
245 if (IS_INF(src))
246 fp_set_nan(dest);
247 /* infinity / anything else = infinity (with approprate sign) */
248 return dest;
249 }
250 if (IS_INF(src)) {
251 /* anything / infinity = zero (with appropriate sign) */
252 dest->exp = 0;
253 dest->mant.m64 = 0;
254 dest->lowmant = 0;
255
256 return dest;
257 }
258
259 /* zeroes */
260 if (IS_ZERO(dest)) {
261 /* zero / zero = NaN */
262 if (IS_ZERO(src))
263 fp_set_nan(dest);
264 /* zero / anything else = zero */
265 return dest;
266 }
267 if (IS_ZERO(src)) {
268 /* anything / zero = infinity (with appropriate sign) */
269 fp_set_sr(FPSR_EXC_DZ);
270 dest->exp = 0x7fff;
271 dest->mant.m64 = 0;
272
273 return dest;
274 }
275
276 exp = dest->exp - src->exp + 0x3fff;
277
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);
285
286 /* now, do the 64-bit divide */
287 fp_dividemant(&temp, dest, src);
288
289 /* normalize it back to 64 bits and stuff it back into the
290 destination struct */
291 if (!temp.m32[0]) {
292 exp--;
293 fp_putmant128(dest, &temp, 32);
294 } else
295 fp_putmant128(dest, &temp, 31);
296
297 if (exp >= 0x7fff) {
298 fp_set_ovrflw(dest);
299 return dest;
300 }
301 dest->exp = exp;
302 if (exp < 0) {
303 fp_set_sr(FPSR_EXC_UNFL);
304 fp_denormalize(dest, -exp);
305 }
306
307 return dest;
308 }
309
310 struct fp_ext *
fp_fsglmul(struct fp_ext * dest,struct fp_ext * src)311 fp_fsglmul(struct fp_ext *dest, struct fp_ext *src)
312 {
313 int exp;
314
315 dprint(PINSTR, "fsglmul\n");
316
317 fp_dyadic_check(dest, src);
318
319 /* calculate the correct sign now, as it's necessary for infinities */
320 dest->sign = src->sign ^ dest->sign;
321
322 /* Handle infinities */
323 if (IS_INF(dest)) {
324 if (IS_ZERO(src))
325 fp_set_nan(dest);
326 return dest;
327 }
328 if (IS_INF(src)) {
329 if (IS_ZERO(dest))
330 fp_set_nan(dest);
331 else
332 fp_copy_ext(dest, src);
333 return dest;
334 }
335
336 /* Of course, as we all know, zero * anything = zero. You may
337 not have known that it might be a positive or negative
338 zero... */
339 if (IS_ZERO(dest) || IS_ZERO(src)) {
340 dest->exp = 0;
341 dest->mant.m64 = 0;
342 dest->lowmant = 0;
343
344 return dest;
345 }
346
347 exp = dest->exp + src->exp - 0x3ffe;
348
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);
353
354 if (exp >= 0x7fff) {
355 fp_set_ovrflw(dest);
356 return dest;
357 }
358 dest->exp = exp;
359 if (exp < 0) {
360 fp_set_sr(FPSR_EXC_UNFL);
361 fp_denormalize(dest, -exp);
362 }
363
364 return dest;
365 }
366
367 struct fp_ext *
fp_fsgldiv(struct fp_ext * dest,struct fp_ext * src)368 fp_fsgldiv(struct fp_ext *dest, struct fp_ext *src)
369 {
370 int exp;
371 unsigned long quot, rem;
372
373 dprint(PINSTR, "fsgldiv\n");
374
375 fp_dyadic_check(dest, src);
376
377 /* calculate the correct sign now, as it's necessary for infinities */
378 dest->sign = src->sign ^ dest->sign;
379
380 /* Handle infinities */
381 if (IS_INF(dest)) {
382 /* infinity / infinity = NaN (quiet, as always) */
383 if (IS_INF(src))
384 fp_set_nan(dest);
385 /* infinity / anything else = infinity (with approprate sign) */
386 return dest;
387 }
388 if (IS_INF(src)) {
389 /* anything / infinity = zero (with appropriate sign) */
390 dest->exp = 0;
391 dest->mant.m64 = 0;
392 dest->lowmant = 0;
393
394 return dest;
395 }
396
397 /* zeroes */
398 if (IS_ZERO(dest)) {
399 /* zero / zero = NaN */
400 if (IS_ZERO(src))
401 fp_set_nan(dest);
402 /* zero / anything else = zero */
403 return dest;
404 }
405 if (IS_ZERO(src)) {
406 /* anything / zero = infinity (with appropriate sign) */
407 fp_set_sr(FPSR_EXC_DZ);
408 dest->exp = 0x7fff;
409 dest->mant.m64 = 0;
410
411 return dest;
412 }
413
414 exp = dest->exp - src->exp + 0x3fff;
415
416 dest->mant.m32[0] &= 0xffffff00;
417 src->mant.m32[0] &= 0xffffff00;
418
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 */
425 } else {
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 */
429 exp--;
430 }
431
432 if (exp >= 0x7fff) {
433 fp_set_ovrflw(dest);
434 return dest;
435 }
436 dest->exp = exp;
437 if (exp < 0) {
438 fp_set_sr(FPSR_EXC_UNFL);
439 fp_denormalize(dest, -exp);
440 }
441
442 return dest;
443 }
444
445 /* fp_roundint: Internal rounding function for use by several of these
446 emulated instructions.
447
448 This one rounds off the fractional part using the rounding mode
449 specified. */
450
fp_roundint(struct fp_ext * dest,int mode)451 static void fp_roundint(struct fp_ext *dest, int mode)
452 {
453 union fp_mant64 oldmant;
454 unsigned long mask;
455
456 if (!fp_normalize_ext(dest))
457 return;
458
459 /* infinities and zeroes */
460 if (IS_INF(dest) || IS_ZERO(dest))
461 return;
462
463 /* first truncate the lower bits */
464 oldmant = dest->mant;
465 switch (dest->exp) {
466 case 0 ... 0x3ffe:
467 dest->mant.m64 = 0;
468 break;
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)
473 return;
474 break;
475 case 0x401f ... 0x403e:
476 dest->mant.m32[1] &= 0xffffffffU << (0x403e - dest->exp);
477 if (oldmant.m32[1] == dest->mant.m32[1])
478 return;
479 break;
480 default:
481 return;
482 }
483 fp_set_sr(FPSR_EXC_INEX2);
484
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.
490
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) */
496
497 switch (mode) {
498 case FPCR_ROUND_RN:
499 switch (dest->exp) {
500 case 0 ... 0x3ffd:
501 return;
502 case 0x3ffe:
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))
508 return;
509 break;
510 case 0x3fff ... 0x401d:
511 mask = 1 << (0x401d - dest->exp);
512 if (!(oldmant.m32[0] & mask))
513 return;
514 if (oldmant.m32[0] & (mask << 1))
515 break;
516 if (!(oldmant.m32[0] << (dest->exp - 0x3ffd)) &&
517 !oldmant.m32[1])
518 return;
519 break;
520 case 0x401e:
521 if (!(oldmant.m32[1] >= 0))
522 return;
523 if (oldmant.m32[0] & 1)
524 break;
525 if (!(oldmant.m32[1] << 1))
526 return;
527 break;
528 case 0x401f ... 0x403d:
529 mask = 1 << (0x403d - dest->exp);
530 if (!(oldmant.m32[1] & mask))
531 return;
532 if (oldmant.m32[1] & (mask << 1))
533 break;
534 if (!(oldmant.m32[1] << (dest->exp - 0x401d)))
535 return;
536 break;
537 default:
538 return;
539 }
540 break;
541 case FPCR_ROUND_RZ:
542 return;
543 default:
544 if (dest->sign ^ (mode - FPCR_ROUND_RM))
545 break;
546 return;
547 }
548
549 switch (dest->exp) {
550 case 0 ... 0x3ffe:
551 dest->exp = 0x3fff;
552 dest->mant.m64 = 1ULL << 63;
553 break;
554 case 0x3fff ... 0x401e:
555 mask = 1 << (0x401e - dest->exp);
556 if (dest->mant.m32[0] += mask)
557 break;
558 dest->mant.m32[0] = 0x80000000;
559 dest->exp++;
560 break;
561 case 0x401f ... 0x403e:
562 mask = 1 << (0x403e - dest->exp);
563 if (dest->mant.m32[1] += mask)
564 break;
565 if (dest->mant.m32[0] += 1)
566 break;
567 dest->mant.m32[0] = 0x80000000;
568 dest->exp++;
569 break;
570 }
571 }
572
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) */
576
577 static struct fp_ext *
modrem_kernel(struct fp_ext * dest,struct fp_ext * src,int mode)578 modrem_kernel(struct fp_ext *dest, struct fp_ext *src, int mode)
579 {
580 struct fp_ext tmp;
581
582 fp_dyadic_check(dest, src);
583
584 /* Infinities and zeros */
585 if (IS_INF(dest) || IS_ZERO(src)) {
586 fp_set_nan(dest);
587 return dest;
588 }
589 if (IS_ZERO(dest) || IS_INF(src))
590 return dest;
591
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);
596 fp_fmul(&tmp, src);
597 fp_fsub(dest, &tmp);
598
599 /* set the quotient byte */
600 fp_set_quotient((dest->mant.m64 & 0x7f) | (dest->sign << 7));
601 return dest;
602 }
603
604 /* fp_fmod: Implements the kernel of the FMOD instruction.
605
606 Again, the argument order is backwards. The result, as defined in
607 the Motorola manuals, is:
608
609 fmod(src,dest) = (dest - (src * floor(dest / src))) */
610
611 struct fp_ext *
fp_fmod(struct fp_ext * dest,struct fp_ext * src)612 fp_fmod(struct fp_ext *dest, struct fp_ext *src)
613 {
614 dprint(PINSTR, "fmod\n");
615 return modrem_kernel(dest, src, FPCR_ROUND_RZ);
616 }
617
618 /* fp_frem: Implements the kernel of the FREM instruction.
619
620 frem(src,dest) = (dest - (src * round(dest / src)))
621 */
622
623 struct fp_ext *
fp_frem(struct fp_ext * dest,struct fp_ext * src)624 fp_frem(struct fp_ext *dest, struct fp_ext *src)
625 {
626 dprint(PINSTR, "frem\n");
627 return modrem_kernel(dest, src, FPCR_ROUND_RN);
628 }
629
630 struct fp_ext *
fp_fint(struct fp_ext * dest,struct fp_ext * src)631 fp_fint(struct fp_ext *dest, struct fp_ext *src)
632 {
633 dprint(PINSTR, "fint\n");
634
635 fp_copy_ext(dest, src);
636
637 fp_roundint(dest, FPDATA->rnd);
638
639 return dest;
640 }
641
642 struct fp_ext *
fp_fintrz(struct fp_ext * dest,struct fp_ext * src)643 fp_fintrz(struct fp_ext *dest, struct fp_ext *src)
644 {
645 dprint(PINSTR, "fintrz\n");
646
647 fp_copy_ext(dest, src);
648
649 fp_roundint(dest, FPCR_ROUND_RZ);
650
651 return dest;
652 }
653
654 struct fp_ext *
fp_fscale(struct fp_ext * dest,struct fp_ext * src)655 fp_fscale(struct fp_ext *dest, struct fp_ext *src)
656 {
657 int scale, oldround;
658
659 dprint(PINSTR, "fscale\n");
660
661 fp_dyadic_check(dest, src);
662
663 /* Infinities */
664 if (IS_INF(src)) {
665 fp_set_nan(dest);
666 return dest;
667 }
668 if (IS_INF(dest))
669 return dest;
670
671 /* zeroes */
672 if (IS_ZERO(src) || IS_ZERO(dest))
673 return dest;
674
675 /* Source exponent out of range */
676 if (src->exp >= 0x400c) {
677 fp_set_ovrflw(dest);
678 return dest;
679 }
680
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;
686
687 /* new exponent */
688 scale += dest->exp;
689
690 if (scale >= 0x7fff) {
691 fp_set_ovrflw(dest);
692 } else if (scale <= 0) {
693 fp_set_sr(FPSR_EXC_UNFL);
694 fp_denormalize(dest, -scale);
695 } else
696 dest->exp = scale;
697
698 return dest;
699 }
700
701