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