[deliverable/linux.git] / arch / mips / math-emu / dp_div.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_div(union ieee754dp x, union ieee754dp y)
{
	u64 rm;
	int re;
	u64 bm;

	COMPXDP;
	COMPYDP;

	EXPLODEXDP;
	EXPLODEYDP;

	ieee754_clearcx();

	FLUSHXDP;
	FLUSHYDP;

	switch (CLPAIR(xc, yc)) {
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
	case CLPAIR(IEEE754_CLASS_SNAN, 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):
	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):
		ieee754_setcx(IEEE754_INVALID_OPERATION);
		return ieee754dp_nanxcpt(ieee754dp_indef());

	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):
		ieee754_setcx(IEEE754_INVALID_OPERATION);
		return ieee754dp_indef();

	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
		return ieee754dp_zero(xs ^ ys);

	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 ieee754dp_inf(xs ^ ys);

	/*
	 * Zero handling
	 */
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
		ieee754_setcx(IEEE754_INVALID_OPERATION);
		return ieee754dp_indef();

	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
		ieee754_setcx(IEEE754_ZERO_DIVIDE);
		return ieee754dp_inf(xs ^ ys);

	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
		return ieee754dp_zero(xs == ys ? 0 : 1);

	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
		DPDNORMX;

	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
		DPDNORMY;
		break;

	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
		DPDNORMX;
		break;

	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
		break;
	}
	assert(xm & DP_HIDDEN_BIT);
	assert(ym & DP_HIDDEN_BIT);

	/* provide rounding space */
	xm <<= 3;
	ym <<= 3;

	/* now the dirty work */

	rm = 0;
	re = xe - ye;

	for (bm = DP_MBIT(DP_FBITS + 2); bm; bm >>= 1) {
		if (xm >= ym) {
			xm -= ym;
			rm |= bm;
			if (xm == 0)
				break;
		}
		xm <<= 1;
	}

	rm <<= 1;
	if (xm)
		rm |= 1;	/* have remainder, set sticky */

	assert(rm);

	/*
	 * Normalise rm to rounding precision ?
	 */
	while ((rm >> (DP_FBITS + 3)) == 0) {
		rm <<= 1;
		re--;
	}

	return ieee754dp_format(xs == ys ? 0 : 1, re, rm);
}
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_div(union ieee754dp x, union ieee754dp y)
1da177e4	25	{
3f7cac41 RB	26	u64 rm;
	27	int re;
	28	u64 bm;
	29
1da177e4 LT	30	COMPXDP;
	31	COMPYDP;
	32
	33	EXPLODEXDP;
	34	EXPLODEYDP;
	35
9e8bad1f	36	ieee754_clearcx();
1da177e4 LT	37
	38	FLUSHXDP;
	39	FLUSHYDP;
	40
	41	switch (CLPAIR(xc, yc)) {
	42	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
	43	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
	44	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
	45	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
	46	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
	47	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
	48	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
	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):
9e8bad1f	53	ieee754_setcx(IEEE754_INVALID_OPERATION);
90efba36	54	return ieee754dp_nanxcpt(ieee754dp_indef());
1da177e4 LT	55
	56	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
	57	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
	58	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
	59	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
	60	return y;
	61
	62	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
	63	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
	64	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
	65	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
	66	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
	67	return x;
	68
	69
3f7cac41 RB	70	/*
	71	* Infinity handling
	72	*/
1da177e4	73	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
9e8bad1f	74	ieee754_setcx(IEEE754_INVALID_OPERATION);
90efba36	75	return ieee754dp_indef();
1da177e4 LT	76
	77	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
	78	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
	79	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
	80	return ieee754dp_zero(xs ^ ys);
	81
	82	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
	83	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
	84	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
	85	return ieee754dp_inf(xs ^ ys);
	86
3f7cac41 RB	87	/*
	88	* Zero handling
	89	*/
1da177e4	90	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
9e8bad1f	91	ieee754_setcx(IEEE754_INVALID_OPERATION);
90efba36	92	return ieee754dp_indef();
1da177e4 LT	93
	94	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
	95	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
9e8bad1f	96	ieee754_setcx(IEEE754_ZERO_DIVIDE);
90efba36	97	return ieee754dp_inf(xs ^ ys);
1da177e4 LT	98
	99	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
	100	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
	101	return ieee754dp_zero(xs == ys ? 0 : 1);
	102
	103	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
	104	DPDNORMX;
	105
	106	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
	107	DPDNORMY;
	108	break;
	109
	110	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
	111	DPDNORMX;
	112	break;
	113
	114	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
	115	break;
	116	}
	117	assert(xm & DP_HIDDEN_BIT);
	118	assert(ym & DP_HIDDEN_BIT);
	119
	120	/* provide rounding space */
	121	xm <<= 3;
	122	ym <<= 3;
	123
3f7cac41	124	/* now the dirty work */
1da177e4	125
3f7cac41 RB	126	rm = 0;
3f7cac41 RB	127	re = xe - ye;
1da177e4	128
3f7cac41 RB	129	for (bm = DP_MBIT(DP_FBITS + 2); bm; bm >>= 1) {
	130	if (xm >= ym) {
	131	xm -= ym;
	132	rm \|= bm;
	133	if (xm == 0)
	134	break;
1da177e4	135	}
3f7cac41 RB	136	xm <<= 1;
	137	}
	138
	139	rm <<= 1;
	140	if (xm)
	141	rm \|= 1; /* have remainder, set sticky */
1da177e4	142
3f7cac41 RB	143	assert(rm);
	144
	145	/*
	146	* Normalise rm to rounding precision ?
	147	*/
	148	while ((rm >> (DP_FBITS + 3)) == 0) {
	149	rm <<= 1;
	150	re--;
1da177e4	151	}
3f7cac41 RB	152
3f7cac41 RB	153	return ieee754dp_format(xs == ys ? 0 : 1, re, rm);
1da177e4	154	}