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