1// SPDX-License-Identifier: GPL-2.0-or-later
2/*
3
4   fp_arith.c: floating-point math routines for the Linux-m68k
5   floating point emulator.
6
7   Copyright (c) 1998-1999 David Huggins-Daines.
8
9   Somewhat based on the AlphaLinux floating point emulator, by David
10   Mosberger-Tang.
11
12 */
13
14#include "fp_emu.h"
15#include "multi_arith.h"
16#include "fp_arith.h"
17
18const struct fp_ext fp_QNaN =
19{
20	.exp = 0x7fff,
21	.mant = { .m64 = ~0 }
22};
23
24const struct fp_ext fp_Inf =
25{
26	.exp = 0x7fff,
27};
28
29/* let's start with the easy ones */
30
31struct fp_ext *fp_fabs(struct fp_ext *dest, struct fp_ext *src)
32{
33	dprint(PINSTR, "fabs\n");
34
35	fp_monadic_check(dest, src);
36
37	dest->sign = 0;
38
39	return dest;
40}
41
42struct fp_ext *fp_fneg(struct fp_ext *dest, struct fp_ext *src)
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, FSADD, FDADD, FSUB,
56   FDSUB, and FCMP instructions. */
57
58struct fp_ext *fp_fadd(struct fp_ext *dest, struct fp_ext *src)
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 != dest->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->sign) {
80				if (FPDATA->rnd == FPCR_ROUND_RM)
81					dest->sign = 1;
82				else
83					dest->sign = 0;
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.m64) {
103			fp_submant(dest, src, dest);
104			dest->sign = !dest->sign;
105		} else
106			fp_submant(dest, dest, src);
107	}
108
109	return dest;
110}
111
112/* fp_fsub: Implements the kernel of the FSUB, FSSUB, and FDSUB
113   instructions.
114
115   Remember that the arguments are in assembler-syntax order! */
116
117struct fp_ext *fp_fsub(struct fp_ext *dest, struct fp_ext *src)
118{
119	dprint(PINSTR, "fsub ");
120
121	src->sign = !src->sign;
122	return fp_fadd(dest, src);
123}
124
125
126struct fp_ext *fp_fcmp(struct fp_ext *dest, struct fp_ext *src)
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
135struct fp_ext *fp_ftst(struct fp_ext *dest, struct fp_ext *src)
136{
137	dprint(PINSTR, "ftst\n");
138
139	(void)dest;
140
141	return src;
142}
143
144struct fp_ext *fp_fmul(struct fp_ext *dest, struct fp_ext *src)
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 it's necessary for infinities */
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 * anything = zero.  You may
171	   not have known that it might be a positive or negative
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 denormalized numbers,
184	   so that the highest bit is set, this makes the
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 expansion */
192	fp_multiplymant(&temp, dest, src);
193
194	/* normalize it back to 64 bits and stuff it back into the
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 FDIV, FSDIV, FDDIV and
216   FSGLDIV instructions.
217
218   Note that the order of the operands is counter-intuitive: instead
219   of src / dest, the result is actually dest / src. */
220
221struct fp_ext *fp_fdiv(struct fp_ext *dest, struct fp_ext *src)
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 it's necessary for infinities */
231	dest->sign = src->sign ^ dest->sign;
232
233	/* Handle infinities */
234	if (IS_INF(dest)) {
235		/* infinity / infinity = NaN (quiet, as always) */
236		if (IS_INF(src))
237			fp_set_nan(dest);
238		/* infinity / anything else = infinity (with appropriate sign) */
239		return dest;
240	}
241	if (IS_INF(src)) {
242		/* anything / infinity = zero (with appropriate sign) */
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 (with appropriate sign) */
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 denormalized numbers,
270	   so that the highest bit is set, this makes lots
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 stuff it back into the
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
301struct fp_ext *fp_fsglmul(struct fp_ext *dest, struct fp_ext *src)
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 it's necessary for infinities */
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 * anything = zero.  You may
327	   not have known that it might be a positive or negative
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.m32[1],
341		 dest->mant.m32[0] & 0xffffff00,
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
357struct fp_ext *fp_fsgldiv(struct fp_ext *dest, struct fp_ext *src)
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 it's necessary for infinities */
367	dest->sign = src->sign ^ dest->sign;
368
369	/* Handle infinities */
370	if (IS_INF(dest)) {
371		/* infinity / infinity = NaN (quiet, as always) */
372		if (IS_INF(src))
373			fp_set_nan(dest);
374		/* infinity / anything else = infinity (with approprate sign) */
375		return dest;
376	}
377	if (IS_INF(src)) {
378		/* anything / infinity = zero (with appropriate sign) */
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 (with appropriate sign) */
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[0]) {
410		fp_sub64(dest->mant, src->mant);
411		fp_div64(quot, rem, dest->mant.m32[0], 0, src->mant.m32[0]);
412		dest->mant.m32[0] = 0x80000000 | (quot >> 1);
413		dest->mant.m32[1] = (quot & 1) | rem;	/* only for rounding */
414	} else {
415		fp_div64(quot, rem, dest->mant.m32[0], 0, src->mant.m32[0]);
416		dest->mant.m32[0] = quot;
417		dest->mant.m32[1] = rem;		/* only for rounding */
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 use by several of these
435   emulated instructions.
436
437   This one rounds off the fractional part using the rounding mode
438   specified. */
439
440static void fp_roundint(struct fp_ext *dest, int mode)
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] &= 0xffffffffU << (0x401e - dest->exp);
460		dest->mant.m32[1] = 0;
461		if (oldmant.m64 == dest->mant.m64)
462			return;
463		break;
464	case 0x401f ... 0x403e:
465		dest->mant.m32[1] &= 0xffffffffU << (0x403e - dest->exp);
466		if (oldmant.m32[1] == dest->mant.m32[1])
467			return;
468		break;
469	default:
470		return;
471	}
472	fp_set_sr(FPSR_EXC_INEX2);
473
474	/* We might want to normalize upwards here... however, since
475	   we know that this is only called on the output of fp_fdiv,
476	   or with the input to fp_fint or fp_fintrz, and the inputs
477	   to all these functions are either normal or denormalized
478	   (no subnormals allowed!), there's really no need.
479
480	   In the case of fp_fdiv, observe that 0x80000000 / 0xffff =
481	   0xffff8000, and the same holds for 128-bit / 64-bit. (i.e. the
482	   smallest possible normal dividend and the largest possible normal
483	   divisor will still produce a normal quotient, therefore, (normal
484	   << 64) / normal is normal in all cases) */
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 input is always normal, so the
493			   guard bit (bit 63) is always set.  therefore, the
494			   only case in which we will NOT round to 1.0 is when
495			   the input is exactly 0.5. */
496			if (oldmant.m64 == (1ULL << 63))
497				return;
498			break;
499		case 0x3fff ... 0x401d:
500			mask = 1 << (0x401d - dest->exp);
501			if (!(oldmant.m32[0] & mask))
502				return;
503			if (oldmant.m32[0] & (mask << 1))
504				break;
505			if (!(oldmant.m32[0] << (dest->exp - 0x3ffd)) &&
506					!oldmant.m32[1])
507				return;
508			break;
509		case 0x401e:
510			if (oldmant.m32[1] & 0x80000000)
511				return;
512			if (oldmant.m32[0] & 1)
513				break;
514			if (!(oldmant.m32[1] << 1))
515				return;
516			break;
517		case 0x401f ... 0x403d:
518			mask = 1 << (0x403d - dest->exp);
519			if (!(oldmant.m32[1] & mask))
520				return;
521			if (oldmant.m32[1] & (mask << 1))
522				break;
523			if (!(oldmant.m32[1] << (dest->exp - 0x401d)))
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_ROUND_RM))
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->exp);
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->exp);
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 and FMOD instructions
563   (which are exactly the same, except for the rounding used on the
564   intermediate value) */
565
566static struct fp_ext *modrem_kernel(struct fp_ext *dest, struct fp_ext *src,
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 smarter way to do this */
582	fp_copy_ext(&tmp, dest);
583	fp_fdiv(&tmp, src);		/* NOTE: src might be modified */
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) | (dest->sign << 7));
590	return dest;
591}
592
593/* fp_fmod: Implements the kernel of the FMOD instruction.
594
595   Again, the argument order is backwards.  The result, as defined in
596   the Motorola manuals, is:
597
598   fmod(src,dest) = (dest - (src * floor(dest / src))) */
599
600struct fp_ext *fp_fmod(struct fp_ext *dest, struct fp_ext *src)
601{
602	dprint(PINSTR, "fmod\n");
603	return modrem_kernel(dest, src, FPCR_ROUND_RZ);
604}
605
606/* fp_frem: Implements the kernel of the FREM instruction.
607
608   frem(src,dest) = (dest - (src * round(dest / src)))
609 */
610
611struct fp_ext *fp_frem(struct fp_ext *dest, struct fp_ext *src)
612{
613	dprint(PINSTR, "frem\n");
614	return modrem_kernel(dest, src, FPCR_ROUND_RN);
615}
616
617struct fp_ext *fp_fint(struct fp_ext *dest, struct fp_ext *src)
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
628struct fp_ext *fp_fintrz(struct fp_ext *dest, struct fp_ext *src)
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
639struct fp_ext *fp_fscale(struct fp_ext *dest, struct fp_ext *src)
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 zero. */
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