[deliverable/linux.git] / arch / mips / math-emu / dp_sub.c

/* IEEE754 floating point arithmetic
 * double precision: common utilities
 */
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.
 *
 *  This program is free software; you can distribute it and/or modify it
 *  under the terms of the GNU General Public License (Version 2) as
 *  published by the Free Software Foundation.
 *
 *  This program is distributed in the hope it will be useful, but WITHOUT
 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 *  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 *  for more details.
 *
 *  You should have received a copy of the GNU General Public License along
 *  with this program; if not, write to the Free Software Foundation, Inc.,
 *  51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA.
 */

#include "ieee754dp.h"

union ieee754dp ieee754dp_sub(union ieee754dp x, union ieee754dp y)
{
	int s;

	COMPXDP;
	COMPYDP;

	EXPLODEXDP;
	EXPLODEYDP;

	ieee754_clearcx();

	FLUSHXDP;
	FLUSHYDP;

	switch (CLPAIR(xc, yc)) {
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
		return ieee754dp_nanxcpt(y);

	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
		return ieee754dp_nanxcpt(x);

	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
		return y;

	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
		return x;


	/*
	 * Infinity handling
	 */
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
		if (xs != ys)
			return x;
		ieee754_setcx(IEEE754_INVALID_OPERATION);
		return ieee754dp_indef();

	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
		return ieee754dp_inf(ys ^ 1);

	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
		return x;

	/*
	 * Zero handling
	 */
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
		if (xs != ys)
			return x;
		else
			return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD);

	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
		return x;

	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
		/* quick fix up */
		DPSIGN(y) ^= 1;
		return y;

	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
		DPDNORMX;
		/* FALL THROUGH */

	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
		/* normalize ym,ye */
		DPDNORMY;
		break;

	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
		/* normalize xm,xe */
		DPDNORMX;
		break;

	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
		break;
	}
	/* flip sign of y and handle as add */
	ys ^= 1;

	assert(xm & DP_HIDDEN_BIT);
	assert(ym & DP_HIDDEN_BIT);


	/* provide guard,round and stick bit dpace */
	xm <<= 3;
	ym <<= 3;

	if (xe > ye) {
		/*
		 * Have to shift y fraction right to align
		 */
		s = xe - ye;
		ym = XDPSRS(ym, s);
		ye += s;
	} else if (ye > xe) {
		/*
		 * Have to shift x fraction right to align
		 */
		s = ye - xe;
		xm = XDPSRS(xm, s);
		xe += s;
	}
	assert(xe == ye);
	assert(xe <= DP_EMAX);

	if (xs == ys) {
		/* generate 28 bit result of adding two 27 bit numbers
		 */
		xm = xm + ym;

		if (xm >> (DP_FBITS + 1 + 3)) { /* carry out */
			xm = XDPSRS1(xm);	/* shift preserving sticky */
			xe++;
		}
	} else {
		if (xm >= ym) {
			xm = xm - ym;
		} else {
			xm = ym - xm;
			xs = ys;
		}
		if (xm == 0) {
			if (ieee754_csr.rm == FPU_CSR_RD)
				return ieee754dp_zero(1);	/* round negative inf. => sign = -1 */
			else
				return ieee754dp_zero(0);	/* other round modes   => sign = 1 */
		}

		/* normalize to rounding precision
		 */
		while ((xm >> (DP_FBITS + 3)) == 0) {
			xm <<= 1;
			xe--;
		}
	}

	return ieee754dp_format(xs, xe, xm);
}
Commit	Line	Data
1da177e4 LT	1	/* IEEE754 floating point arithmetic
	2	* double precision: common utilities
	3	*/
	4	/*
	5	* MIPS floating point support
	6	* Copyright (C) 1994-2000 Algorithmics Ltd.
1da177e4	7	*
1da177e4 LT	8	* This program is free software; you can distribute it and/or modify it
	9	* under the terms of the GNU General Public License (Version 2) as
	10	* published by the Free Software Foundation.
	11	*
	12	* This program is distributed in the hope it will be useful, but WITHOUT
	13	* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
	14	* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
	15	* for more details.
	16	*
	17	* You should have received a copy of the GNU General Public License along
	18	* with this program; if not, write to the Free Software Foundation, Inc.,
3f7cac41	19	* 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
1da177e4 LT	20	*/
1da177e4 LT	21
1da177e4 LT	22	#include "ieee754dp.h"
1da177e4 LT	23
2209bcb1	24	union ieee754dp ieee754dp_sub(union ieee754dp x, union ieee754dp y)
1da177e4	25	{
3f7cac41 RB	26	int s;
3f7cac41 RB	27
1da177e4 LT	28	COMPXDP;
	29	COMPYDP;
	30
	31	EXPLODEXDP;
	32	EXPLODEYDP;
	33
9e8bad1f	34	ieee754_clearcx();
1da177e4 LT	35
	36	FLUSHXDP;
	37	FLUSHYDP;
	38
	39	switch (CLPAIR(xc, yc)) {
1da177e4	40	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
1da177e4 LT	41	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
	42	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
	43	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
	44	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
d5afa7e9 MR	45	return ieee754dp_nanxcpt(y);
	46
	47	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
	48	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
1da177e4 LT	49	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
	50	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
	51	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
	52	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
d5afa7e9	53	return ieee754dp_nanxcpt(x);
1da177e4 LT	54
	55	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
	56	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
	57	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
	58	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
	59	return y;
	60
	61	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
	62	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
	63	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
	64	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
	65	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
	66	return x;
	67
	68
3f7cac41 RB	69	/*
	70	* Infinity handling
	71	*/
1da177e4 LT	72	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
	73	if (xs != ys)
	74	return x;
9e8bad1f	75	ieee754_setcx(IEEE754_INVALID_OPERATION);
90efba36	76	return ieee754dp_indef();
1da177e4 LT	77
	78	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
	79	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
	80	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
	81	return ieee754dp_inf(ys ^ 1);
	82
	83	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
	84	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
	85	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
	86	return x;
	87
3f7cac41 RB	88	/*
	89	* Zero handling
	90	*/
1da177e4 LT	91	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
	92	if (xs != ys)
	93	return x;
	94	else
56a64733	95	return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD);
1da177e4 LT	96
	97	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
	98	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
	99	return x;
	100
	101	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
	102	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
	103	/* quick fix up */
	104	DPSIGN(y) ^= 1;
	105	return y;
	106
	107	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
	108	DPDNORMX;
af901ca1	109	/* FALL THROUGH */
1da177e4 LT	110
	111	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
	112	/* normalize ym,ye */
	113	DPDNORMY;
	114	break;
	115
	116	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
	117	/* normalize xm,xe */
	118	DPDNORMX;
	119	break;
	120
	121	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
	122	break;
	123	}
	124	/* flip sign of y and handle as add */
	125	ys ^= 1;
	126
	127	assert(xm & DP_HIDDEN_BIT);
	128	assert(ym & DP_HIDDEN_BIT);
	129
	130
	131	/* provide guard,round and stick bit dpace */
	132	xm <<= 3;
	133	ym <<= 3;
	134
	135	if (xe > ye) {
3f7cac41 RB	136	/*
3f7cac41 RB	137	* Have to shift y fraction right to align
1da177e4	138	*/
3f7cac41	139	s = xe - ye;
1da177e4 LT	140	ym = XDPSRS(ym, s);
	141	ye += s;
	142	} else if (ye > xe) {
3f7cac41 RB	143	/*
3f7cac41 RB	144	* Have to shift x fraction right to align
1da177e4	145	*/
3f7cac41	146	s = ye - xe;
1da177e4 LT	147	xm = XDPSRS(xm, s);
	148	xe += s;
	149	}
	150	assert(xe == ye);
	151	assert(xe <= DP_EMAX);
	152
	153	if (xs == ys) {
	154	/* generate 28 bit result of adding two 27 bit numbers
	155	*/
	156	xm = xm + ym;
1da177e4	157
ad8fb553	158	if (xm >> (DP_FBITS + 1 + 3)) { /* carry out */
1da177e4 LT	159	xm = XDPSRS1(xm); /* shift preserving sticky */
	160	xe++;
	161	}
	162	} else {
	163	if (xm >= ym) {
	164	xm = xm - ym;
1da177e4 LT	165	} else {
1da177e4 LT	166	xm = ym - xm;
1da177e4 LT	167	xs = ys;
	168	}
	169	if (xm == 0) {
56a64733	170	if (ieee754_csr.rm == FPU_CSR_RD)
1da177e4 LT	171	return ieee754dp_zero(1); /* round negative inf. => sign = -1 */
	172	else
	173	return ieee754dp_zero(0); /* other round modes => sign = 1 */
	174	}
	175
	176	/* normalize to rounding precision
	177	*/
ad8fb553	178	while ((xm >> (DP_FBITS + 3)) == 0) {
1da177e4 LT	179	xm <<= 1;
	180	xe--;
	181	}
	182	}
90efba36 RB	183
90efba36 RB	184	return ieee754dp_format(xs, xe, xm);
1da177e4	185	}