[linux-2.6-block.git] / arch / m68k / math-emu / fp_arith.c

// SPDX-License-Identifier: GPL-2.0-or-later
/*

   fp_arith.c: floating-point math routines for the Linux-m68k
   floating point emulator.

   Copyright (c) 1998-1999 David Huggins-Daines.

   Somewhat based on the AlphaLinux floating point emulator, by David
   Mosberger-Tang.

 */

#include "fp_emu.h"
#include "multi_arith.h"
#include "fp_arith.h"

const struct fp_ext fp_QNaN =
{
	.exp = 0x7fff,
	.mant = { .m64 = ~0 }
};

const struct fp_ext fp_Inf =
{
	.exp = 0x7fff,
};

/* let's start with the easy ones */

struct fp_ext *
fp_fabs(struct fp_ext *dest, struct fp_ext *src)
{
	dprint(PINSTR, "fabs\n");

	fp_monadic_check(dest, src);

	dest->sign = 0;

	return dest;
}

struct fp_ext *
fp_fneg(struct fp_ext *dest, struct fp_ext *src)
{
	dprint(PINSTR, "fneg\n");

	fp_monadic_check(dest, src);

	dest->sign = !dest->sign;

	return dest;
}

/* Now, the slightly harder ones */

/* fp_fadd: Implements the kernel of the FADD, FSADD, FDADD, FSUB,
   FDSUB, and FCMP instructions. */

struct fp_ext *
fp_fadd(struct fp_ext *dest, struct fp_ext *src)
{
	int diff;

	dprint(PINSTR, "fadd\n");

	fp_dyadic_check(dest, src);

	if (IS_INF(dest)) {
		/* infinity - infinity == NaN */
		if (IS_INF(src) && (src->sign != dest->sign))
			fp_set_nan(dest);
		return dest;
	}
	if (IS_INF(src)) {
		fp_copy_ext(dest, src);
		return dest;
	}

	if (IS_ZERO(dest)) {
		if (IS_ZERO(src)) {
			if (src->sign != dest->sign) {
				if (FPDATA->rnd == FPCR_ROUND_RM)
					dest->sign = 1;
				else
					dest->sign = 0;
			}
		} else
			fp_copy_ext(dest, src);
		return dest;
	}

	dest->lowmant = src->lowmant = 0;

	if ((diff = dest->exp - src->exp) > 0)
		fp_denormalize(src, diff);
	else if ((diff = -diff) > 0)
		fp_denormalize(dest, diff);

	if (dest->sign == src->sign) {
		if (fp_addmant(dest, src))
			if (!fp_addcarry(dest))
				return dest;
	} else {
		if (dest->mant.m64 < src->mant.m64) {
			fp_submant(dest, src, dest);
			dest->sign = !dest->sign;
		} else
			fp_submant(dest, dest, src);
	}

	return dest;
}

/* fp_fsub: Implements the kernel of the FSUB, FSSUB, and FDSUB
   instructions.

   Remember that the arguments are in assembler-syntax order! */

struct fp_ext *
fp_fsub(struct fp_ext *dest, struct fp_ext *src)
{
	dprint(PINSTR, "fsub ");

	src->sign = !src->sign;
	return fp_fadd(dest, src);
}


struct fp_ext *
fp_fcmp(struct fp_ext *dest, struct fp_ext *src)
{
	dprint(PINSTR, "fcmp ");

	FPDATA->temp[1] = *dest;
	src->sign = !src->sign;
	return fp_fadd(&FPDATA->temp[1], src);
}

struct fp_ext *
fp_ftst(struct fp_ext *dest, struct fp_ext *src)
{
	dprint(PINSTR, "ftst\n");

	(void)dest;

	return src;
}

struct fp_ext *
fp_fmul(struct fp_ext *dest, struct fp_ext *src)
{
	union fp_mant128 temp;
	int exp;

	dprint(PINSTR, "fmul\n");

	fp_dyadic_check(dest, src);

	/* calculate the correct sign now, as it's necessary for infinities */
	dest->sign = src->sign ^ dest->sign;

	/* Handle infinities */
	if (IS_INF(dest)) {
		if (IS_ZERO(src))
			fp_set_nan(dest);
		return dest;
	}
	if (IS_INF(src)) {
		if (IS_ZERO(dest))
			fp_set_nan(dest);
		else
			fp_copy_ext(dest, src);
		return dest;
	}

	/* Of course, as we all know, zero * anything = zero.  You may
	   not have known that it might be a positive or negative
	   zero... */
	if (IS_ZERO(dest) || IS_ZERO(src)) {
		dest->exp = 0;
		dest->mant.m64 = 0;
		dest->lowmant = 0;

		return dest;
	}

	exp = dest->exp + src->exp - 0x3ffe;

	/* shift up the mantissa for denormalized numbers,
	   so that the highest bit is set, this makes the
	   shift of the result below easier */
	if ((long)dest->mant.m32[0] >= 0)
		exp -= fp_overnormalize(dest);
	if ((long)src->mant.m32[0] >= 0)
		exp -= fp_overnormalize(src);

	/* now, do a 64-bit multiply with expansion */
	fp_multiplymant(&temp, dest, src);

	/* normalize it back to 64 bits and stuff it back into the
	   destination struct */
	if ((long)temp.m32[0] > 0) {
		exp--;
		fp_putmant128(dest, &temp, 1);
	} else
		fp_putmant128(dest, &temp, 0);

	if (exp >= 0x7fff) {
		fp_set_ovrflw(dest);
		return dest;
	}
	dest->exp = exp;
	if (exp < 0) {
		fp_set_sr(FPSR_EXC_UNFL);
		fp_denormalize(dest, -exp);
	}

	return dest;
}

/* fp_fdiv: Implements the "kernel" of the FDIV, FSDIV, FDDIV and
   FSGLDIV instructions.

   Note that the order of the operands is counter-intuitive: instead
   of src / dest, the result is actually dest / src. */

struct fp_ext *
fp_fdiv(struct fp_ext *dest, struct fp_ext *src)
{
	union fp_mant128 temp;
	int exp;

	dprint(PINSTR, "fdiv\n");

	fp_dyadic_check(dest, src);

	/* calculate the correct sign now, as it's necessary for infinities */
	dest->sign = src->sign ^ dest->sign;

	/* Handle infinities */
	if (IS_INF(dest)) {
		/* infinity / infinity = NaN (quiet, as always) */
		if (IS_INF(src))
			fp_set_nan(dest);
		/* infinity / anything else = infinity (with approprate sign) */
		return dest;
	}
	if (IS_INF(src)) {
		/* anything / infinity = zero (with appropriate sign) */
		dest->exp = 0;
		dest->mant.m64 = 0;
		dest->lowmant = 0;

		return dest;
	}

	/* zeroes */
	if (IS_ZERO(dest)) {
		/* zero / zero = NaN */
		if (IS_ZERO(src))
			fp_set_nan(dest);
		/* zero / anything else = zero */
		return dest;
	}
	if (IS_ZERO(src)) {
		/* anything / zero = infinity (with appropriate sign) */
		fp_set_sr(FPSR_EXC_DZ);
		dest->exp = 0x7fff;
		dest->mant.m64 = 0;

		return dest;
	}

	exp = dest->exp - src->exp + 0x3fff;

	/* shift up the mantissa for denormalized numbers,
	   so that the highest bit is set, this makes lots
	   of things below easier */
	if ((long)dest->mant.m32[0] >= 0)
		exp -= fp_overnormalize(dest);
	if ((long)src->mant.m32[0] >= 0)
		exp -= fp_overnormalize(src);

	/* now, do the 64-bit divide */
	fp_dividemant(&temp, dest, src);

	/* normalize it back to 64 bits and stuff it back into the
	   destination struct */
	if (!temp.m32[0]) {
		exp--;
		fp_putmant128(dest, &temp, 32);
	} else
		fp_putmant128(dest, &temp, 31);

	if (exp >= 0x7fff) {
		fp_set_ovrflw(dest);
		return dest;
	}
	dest->exp = exp;
	if (exp < 0) {
		fp_set_sr(FPSR_EXC_UNFL);
		fp_denormalize(dest, -exp);
	}

	return dest;
}

struct fp_ext *
fp_fsglmul(struct fp_ext *dest, struct fp_ext *src)
{
	int exp;

	dprint(PINSTR, "fsglmul\n");

	fp_dyadic_check(dest, src);

	/* calculate the correct sign now, as it's necessary for infinities */
	dest->sign = src->sign ^ dest->sign;

	/* Handle infinities */
	if (IS_INF(dest)) {
		if (IS_ZERO(src))
			fp_set_nan(dest);
		return dest;
	}
	if (IS_INF(src)) {
		if (IS_ZERO(dest))
			fp_set_nan(dest);
		else
			fp_copy_ext(dest, src);
		return dest;
	}

	/* Of course, as we all know, zero * anything = zero.  You may
	   not have known that it might be a positive or negative
	   zero... */
	if (IS_ZERO(dest) || IS_ZERO(src)) {
		dest->exp = 0;
		dest->mant.m64 = 0;
		dest->lowmant = 0;

		return dest;
	}

	exp = dest->exp + src->exp - 0x3ffe;

	/* do a 32-bit multiply */
	fp_mul64(dest->mant.m32[0], dest->mant.m32[1],
		 dest->mant.m32[0] & 0xffffff00,
		 src->mant.m32[0] & 0xffffff00);

	if (exp >= 0x7fff) {
		fp_set_ovrflw(dest);
		return dest;
	}
	dest->exp = exp;
	if (exp < 0) {
		fp_set_sr(FPSR_EXC_UNFL);
		fp_denormalize(dest, -exp);
	}

	return dest;
}

struct fp_ext *
fp_fsgldiv(struct fp_ext *dest, struct fp_ext *src)
{
	int exp;
	unsigned long quot, rem;

	dprint(PINSTR, "fsgldiv\n");

	fp_dyadic_check(dest, src);

	/* calculate the correct sign now, as it's necessary for infinities */
	dest->sign = src->sign ^ dest->sign;

	/* Handle infinities */
	if (IS_INF(dest)) {
		/* infinity / infinity = NaN (quiet, as always) */
		if (IS_INF(src))
			fp_set_nan(dest);
		/* infinity / anything else = infinity (with approprate sign) */
		return dest;
	}
	if (IS_INF(src)) {
		/* anything / infinity = zero (with appropriate sign) */
		dest->exp = 0;
		dest->mant.m64 = 0;
		dest->lowmant = 0;

		return dest;
	}

	/* zeroes */
	if (IS_ZERO(dest)) {
		/* zero / zero = NaN */
		if (IS_ZERO(src))
			fp_set_nan(dest);
		/* zero / anything else = zero */
		return dest;
	}
	if (IS_ZERO(src)) {
		/* anything / zero = infinity (with appropriate sign) */
		fp_set_sr(FPSR_EXC_DZ);
		dest->exp = 0x7fff;
		dest->mant.m64 = 0;

		return dest;
	}

	exp = dest->exp - src->exp + 0x3fff;

	dest->mant.m32[0] &= 0xffffff00;
	src->mant.m32[0] &= 0xffffff00;

	/* do the 32-bit divide */
	if (dest->mant.m32[0] >= src->mant.m32[0]) {
		fp_sub64(dest->mant, src->mant);
		fp_div64(quot, rem, dest->mant.m32[0], 0, src->mant.m32[0]);
		dest->mant.m32[0] = 0x80000000 | (quot >> 1);
		dest->mant.m32[1] = (quot & 1) | rem;	/* only for rounding */
	} else {
		fp_div64(quot, rem, dest->mant.m32[0], 0, src->mant.m32[0]);
		dest->mant.m32[0] = quot;
		dest->mant.m32[1] = rem;		/* only for rounding */
		exp--;
	}

	if (exp >= 0x7fff) {
		fp_set_ovrflw(dest);
		return dest;
	}
	dest->exp = exp;
	if (exp < 0) {
		fp_set_sr(FPSR_EXC_UNFL);
		fp_denormalize(dest, -exp);
	}

	return dest;
}

/* fp_roundint: Internal rounding function for use by several of these
   emulated instructions.

   This one rounds off the fractional part using the rounding mode
   specified. */

static void fp_roundint(struct fp_ext *dest, int mode)
{
	union fp_mant64 oldmant;
	unsigned long mask;

	if (!fp_normalize_ext(dest))
		return;

	/* infinities and zeroes */
	if (IS_INF(dest) || IS_ZERO(dest))
		return;

	/* first truncate the lower bits */
	oldmant = dest->mant;
	switch (dest->exp) {
	case 0 ... 0x3ffe:
		dest->mant.m64 = 0;
		break;
	case 0x3fff ... 0x401e:
		dest->mant.m32[0] &= 0xffffffffU << (0x401e - dest->exp);
		dest->mant.m32[1] = 0;
		if (oldmant.m64 == dest->mant.m64)
			return;
		break;
	case 0x401f ... 0x403e:
		dest->mant.m32[1] &= 0xffffffffU << (0x403e - dest->exp);
		if (oldmant.m32[1] == dest->mant.m32[1])
			return;
		break;
	default:
		return;
	}
	fp_set_sr(FPSR_EXC_INEX2);

	/* We might want to normalize upwards here... however, since
	   we know that this is only called on the output of fp_fdiv,
	   or with the input to fp_fint or fp_fintrz, and the inputs
	   to all these functions are either normal or denormalized
	   (no subnormals allowed!), there's really no need.

	   In the case of fp_fdiv, observe that 0x80000000 / 0xffff =
	   0xffff8000, and the same holds for 128-bit / 64-bit. (i.e. the
	   smallest possible normal dividend and the largest possible normal
	   divisor will still produce a normal quotient, therefore, (normal
	   << 64) / normal is normal in all cases) */

	switch (mode) {
	case FPCR_ROUND_RN:
		switch (dest->exp) {
		case 0 ... 0x3ffd:
			return;
		case 0x3ffe:
			/* As noted above, the input is always normal, so the
			   guard bit (bit 63) is always set.  therefore, the
			   only case in which we will NOT round to 1.0 is when
			   the input is exactly 0.5. */
			if (oldmant.m64 == (1ULL << 63))
				return;
			break;
		case 0x3fff ... 0x401d:
			mask = 1 << (0x401d - dest->exp);
			if (!(oldmant.m32[0] & mask))
				return;
			if (oldmant.m32[0] & (mask << 1))
				break;
			if (!(oldmant.m32[0] << (dest->exp - 0x3ffd)) &&
					!oldmant.m32[1])
				return;
			break;
		case 0x401e:
			if (oldmant.m32[1] & 0x80000000)
				return;
			if (oldmant.m32[0] & 1)
				break;
			if (!(oldmant.m32[1] << 1))
				return;
			break;
		case 0x401f ... 0x403d:
			mask = 1 << (0x403d - dest->exp);
			if (!(oldmant.m32[1] & mask))
				return;
			if (oldmant.m32[1] & (mask << 1))
				break;
			if (!(oldmant.m32[1] << (dest->exp - 0x401d)))
				return;
			break;
		default:
			return;
		}
		break;
	case FPCR_ROUND_RZ:
		return;
	default:
		if (dest->sign ^ (mode - FPCR_ROUND_RM))
			break;
		return;
	}

	switch (dest->exp) {
	case 0 ... 0x3ffe:
		dest->exp = 0x3fff;
		dest->mant.m64 = 1ULL << 63;
		break;
	case 0x3fff ... 0x401e:
		mask = 1 << (0x401e - dest->exp);
		if (dest->mant.m32[0] += mask)
			break;
		dest->mant.m32[0] = 0x80000000;
		dest->exp++;
		break;
	case 0x401f ... 0x403e:
		mask = 1 << (0x403e - dest->exp);
		if (dest->mant.m32[1] += mask)
			break;
		if (dest->mant.m32[0] += 1)
                        break;
		dest->mant.m32[0] = 0x80000000;
                dest->exp++;
		break;
	}
}

/* modrem_kernel: Implementation of the FREM and FMOD instructions
   (which are exactly the same, except for the rounding used on the
   intermediate value) */

static struct fp_ext *
modrem_kernel(struct fp_ext *dest, struct fp_ext *src, int mode)
{
	struct fp_ext tmp;

	fp_dyadic_check(dest, src);

	/* Infinities and zeros */
	if (IS_INF(dest) || IS_ZERO(src)) {
		fp_set_nan(dest);
		return dest;
	}
	if (IS_ZERO(dest) || IS_INF(src))
		return dest;

	/* FIXME: there is almost certainly a smarter way to do this */
	fp_copy_ext(&tmp, dest);
	fp_fdiv(&tmp, src);		/* NOTE: src might be modified */
	fp_roundint(&tmp, mode);
	fp_fmul(&tmp, src);
	fp_fsub(dest, &tmp);

	/* set the quotient byte */
	fp_set_quotient((dest->mant.m64 & 0x7f) | (dest->sign << 7));
	return dest;
}

/* fp_fmod: Implements the kernel of the FMOD instruction.

   Again, the argument order is backwards.  The result, as defined in
   the Motorola manuals, is:

   fmod(src,dest) = (dest - (src * floor(dest / src))) */

struct fp_ext *
fp_fmod(struct fp_ext *dest, struct fp_ext *src)
{
	dprint(PINSTR, "fmod\n");
	return modrem_kernel(dest, src, FPCR_ROUND_RZ);
}

/* fp_frem: Implements the kernel of the FREM instruction.

   frem(src,dest) = (dest - (src * round(dest / src)))
 */

struct fp_ext *
fp_frem(struct fp_ext *dest, struct fp_ext *src)
{
	dprint(PINSTR, "frem\n");
	return modrem_kernel(dest, src, FPCR_ROUND_RN);
}

struct fp_ext *
fp_fint(struct fp_ext *dest, struct fp_ext *src)
{
	dprint(PINSTR, "fint\n");

	fp_copy_ext(dest, src);

	fp_roundint(dest, FPDATA->rnd);

	return dest;
}

struct fp_ext *
fp_fintrz(struct fp_ext *dest, struct fp_ext *src)
{
	dprint(PINSTR, "fintrz\n");

	fp_copy_ext(dest, src);

	fp_roundint(dest, FPCR_ROUND_RZ);

	return dest;
}

struct fp_ext *
fp_fscale(struct fp_ext *dest, struct fp_ext *src)
{
	int scale, oldround;

	dprint(PINSTR, "fscale\n");

	fp_dyadic_check(dest, src);

	/* Infinities */
	if (IS_INF(src)) {
		fp_set_nan(dest);
		return dest;
	}
	if (IS_INF(dest))
		return dest;

	/* zeroes */
	if (IS_ZERO(src) || IS_ZERO(dest))
		return dest;

	/* Source exponent out of range */
	if (src->exp >= 0x400c) {
		fp_set_ovrflw(dest);
		return dest;
	}

	/* src must be rounded with round to zero. */
	oldround = FPDATA->rnd;
	FPDATA->rnd = FPCR_ROUND_RZ;
	scale = fp_conv_ext2long(src);
	FPDATA->rnd = oldround;

	/* new exponent */
	scale += dest->exp;

	if (scale >= 0x7fff) {
		fp_set_ovrflw(dest);
	} else if (scale <= 0) {
		fp_set_sr(FPSR_EXC_UNFL);
		fp_denormalize(dest, -scale);
	} else
		dest->exp = scale;

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