[Linux-mips-commits]CVS: linux/arch/mips/math-emu dp_add.c,NONE,1.1 dp_cmp.c,NONE,1.1 dp_div.c,NONE,1.1 dp_fsp.c,NONE,1.1 dp_mul.c,NONE,1.1 dp_simple.c,NONE,1.1 dp_sqrt.c,NONE,1.1 dp_sub.c,1.2,1.3 ieee754.c,NONE,1.1 ieee754.h,NONE,1.1 ieee754dp.c,NONE,1.1 ieee754int.h,NONE,1.1 ieee754sp.c,NONE,1.1 kernel_linkage.c,NONE,1.1 sp_add.c,NONE,1.1 sp_cmp.c,NONE,1.1 sp_div.c,NONE,1.1 sp_fdp.c,NONE,1.1 sp_mul.c,NONE,1.1 sp_simple.c,NONE,1.1 sp_sqrt.c,NONE,1.1 sp_sub.c,NONE,1.1 cp1emu.c,1.6,1.7

[James S. <jsi...@us...>]

Update of /cvsroot/linux-mips/linux/arch/mips/math-emu
In directory usw-pr-cvs1:/tmp/cvs-serv23991

Modified Files:
	cp1emu.c 
Added Files:
	dp_add.c dp_cmp.c dp_div.c dp_fsp.c dp_mul.c dp_simple.c 
	dp_sqrt.c dp_sub.c ieee754.c ieee754.h ieee754dp.c 
	ieee754int.h ieee754sp.c kernel_linkage.c sp_add.c sp_cmp.c 
	sp_div.c sp_fdp.c sp_mul.c sp_simple.c sp_sqrt.c sp_sub.c 
Log Message:
Lots of FPU bug fixes from Kjeld Borch Egevang.

--- NEW FILE: dp_add.c ---
/* IEEE754 floating point arithmetic
 * double precision: common utilities
 */
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.  All rights reserved.
 * http://www.algor.co.uk
 *
 * ########################################################################
 *
 *  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.,
 *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
 *
 * ########################################################################
 *
 */


#include "ieee754dp.h"

ieee754dp ieee754dp_add(ieee754dp x, ieee754dp y)
{
	COMPXDP;
	COMPYDP;

	EXPLODEXDP;
	EXPLODEYDP;

	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):
		SETCX(IEEE754_INVALID_OPERATION);
		return ieee754dp_nanxcpt(ieee754dp_indef(), "add", x, y);

	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;


		/* Inifity handeling 
		 */

	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
		if (xs == ys)
			return x;
		SETCX(IEEE754_INVALID_OPERATION);
		return ieee754dp_xcpt(ieee754dp_indef(), "add", x, y);

	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 y;

	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 handeling 
		 */

	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
		if (xs == ys)
			return x;
		else
			return ieee754dp_zero(ieee754_csr.rm ==
					      IEEE754_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):
		return y;

	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
		DPDNORMX;

		/* FALL THROUGH */

	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 guard,round and stick bit space */
	xm <<= 3;
	ym <<= 3;

	if (xe > ye) {
		/* have to shift y fraction right to align
		 */
		int s = xe - ye;
		ym = XDPSRS(ym, s);
		ye += s;
	} else if (ye > xe) {
		/* have to shift x fraction right to align
		 */
		int 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
		 * leaving result in xm,xs,xe
		 */
		xm = xm + ym;
		xe = xe;
		xs = xs;

		if (xm >> (DP_MBITS + 1 + 3)) {	/* carry out */
			xm = XDPSRS1(xm);
			xe++;
		}
	} else {
		if (xm >= ym) {
			xm = xm - ym;
			xe = xe;
			xs = xs;
		} else {
			xm = ym - xm;
			xe = xe;
			xs = ys;
		}
		if (xm == 0)
			return ieee754dp_zero(ieee754_csr.rm ==
					      IEEE754_RD);

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

	}
	DPNORMRET2(xs, xe, xm, "add", x, y);
}

--- NEW FILE: dp_cmp.c ---
/* IEEE754 floating point arithmetic
 * double precision: common utilities
 */
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.  All rights reserved.
 * http://www.algor.co.uk
 *
 * ########################################################################
 *
 *  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.,
 *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
 *
 * ########################################################################
 */


#include "ieee754dp.h"

int ieee754dp_cmp(ieee754dp x, ieee754dp y, int cmp)
{
	COMPXDP;
	COMPYDP;

	EXPLODEXDP;
	EXPLODEYDP;
	FLUSHXDP;
	FLUSHYDP;
	CLEARCX;	/* Even clear inexact flag here */

	if (ieee754dp_isnan(x) || ieee754dp_isnan(y)) {
		if (xc == IEEE754_CLASS_SNAN || yc == IEEE754_CLASS_SNAN)
			SETCX(IEEE754_INVALID_OPERATION);
		if (cmp & IEEE754_CUN)
			return 1;
		if (cmp & (IEEE754_CLT | IEEE754_CGT)) {
			if (SETCX(IEEE754_INVALID_OPERATION))
				return ieee754si_xcpt(0, "fcmpf", x);
		}
		return 0;
	} else {
		long long int vx = x.bits;
		long long int vy = y.bits;

		if (vx < 0)
			vx = -vx ^ DP_SIGN_BIT;
		if (vy < 0)
			vy = -vy ^ DP_SIGN_BIT;

		if (vx < vy)
			return (cmp & IEEE754_CLT) != 0;
		else if (vx == vy)
			return (cmp & IEEE754_CEQ) != 0;
		else
			return (cmp & IEEE754_CGT) != 0;
	}
}

--- NEW FILE: dp_div.c ---
/* IEEE754 floating point arithmetic
 * double precision: common utilities
 */
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.  All rights reserved.
 * http://www.algor.co.uk
 *
 * ########################################################################
 *
 *  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.,
 *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
 *
 * ########################################################################
 */


#include "ieee754dp.h"

ieee754dp ieee754dp_div(ieee754dp x, ieee754dp y)
{
	COMPXDP;
	COMPYDP;

	EXPLODEXDP;
	EXPLODEYDP;

	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):
		SETCX(IEEE754_INVALID_OPERATION);
		return ieee754dp_nanxcpt(ieee754dp_indef(), "div", x, y);

	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 handeling 
		 */

	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
		SETCX(IEEE754_INVALID_OPERATION);
		return ieee754dp_xcpt(ieee754dp_indef(), "div", x, y);

	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 handeling 
		 */

	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
		SETCX(IEEE754_INVALID_OPERATION);
		return ieee754dp_xcpt(ieee754dp_indef(), "div", x, y);

	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
		SETCX(IEEE754_ZERO_DIVIDE);
		return ieee754dp_xcpt(ieee754dp_inf(xs ^ ys), "div", x, y);

	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 */

		unsigned long long rm = 0;
		int re = xe - ye;
		unsigned long long bm;

		for (bm = DP_MBIT(DP_MBITS + 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_MBITS + 3)) == 0) {
			rm <<= 1;
			re--;
		}

		DPNORMRET2(xs == ys ? 0 : 1, re, rm, "div", x, y);
	}
}

--- NEW FILE: dp_fsp.c ---
/* IEEE754 floating point arithmetic
 * double precision: common utilities
 */
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.  All rights reserved.
 * http://www.algor.co.uk
 *
 * ########################################################################
 *
 *  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.,
 *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
 *
 * ########################################################################
 */


#include "ieee754dp.h"

ieee754dp ieee754dp_fsp(ieee754sp x)
{
	COMPXSP;

	EXPLODEXSP;

	CLEARCX;

	FLUSHXSP;

	switch (xc) {
	case IEEE754_CLASS_SNAN:
		SETCX(IEEE754_INVALID_OPERATION);
		return ieee754dp_nanxcpt(ieee754dp_indef(), "fsp");
	case IEEE754_CLASS_QNAN:
		return ieee754dp_nanxcpt(builddp(xs,
						 DP_EMAX + 1 + DP_EBIAS,
						 ((unsigned long long) xm
						  << (DP_MBITS -
						      SP_MBITS))), "fsp",
					 x);
	case IEEE754_CLASS_INF:
		return ieee754dp_inf(xs);
	case IEEE754_CLASS_ZERO:
		return ieee754dp_zero(xs);
	case IEEE754_CLASS_DNORM:
		/* normalize */
		while ((xm >> SP_MBITS) == 0) {
			xm <<= 1;
			xe--;
		}
		break;
	case IEEE754_CLASS_NORM:
		break;
	}

	/* CANT possibly overflow,underflow, or need rounding
	 */

	/* drop the hidden bit */
	xm &= ~SP_HIDDEN_BIT;

	return builddp(xs, xe + DP_EBIAS,
		       (unsigned long long) xm << (DP_MBITS - SP_MBITS));
}

--- NEW FILE: dp_mul.c ---
/* IEEE754 floating point arithmetic
 * double precision: common utilities
 */
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.  All rights reserved.
 * http://www.algor.co.uk
 *
 * ########################################################################
 *
 *  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.,
 *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
 *
 * ########################################################################
 */


#include "ieee754dp.h"

ieee754dp ieee754dp_mul(ieee754dp x, ieee754dp y)
{
	COMPXDP;
	COMPYDP;

	EXPLODEXDP;
	EXPLODEYDP;

	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):
		SETCX(IEEE754_INVALID_OPERATION);
		return ieee754dp_nanxcpt(ieee754dp_indef(), "mul", x, y);

	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 handeling */

	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
		SETCX(IEEE754_INVALID_OPERATION);
		return ieee754dp_xcpt(ieee754dp_indef(), "mul", x, y);

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

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


	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;
	}
	/* rm = xm * ym, re = xe+ye basicly */
	assert(xm & DP_HIDDEN_BIT);
	assert(ym & DP_HIDDEN_BIT);
	{
		int re = xe + ye;
		int rs = xs ^ ys;
		unsigned long long rm;

		/* shunt to top of word */
		xm <<= 64 - (DP_MBITS + 1);
		ym <<= 64 - (DP_MBITS + 1);

		/* multiply 32bits xm,ym to give high 32bits rm with stickness
		 */

		/* 32 * 32 => 64 */
#define DPXMULT(x,y)	((unsigned long long)(x) * (unsigned long long)y)

		{
			unsigned lxm = xm;
			unsigned hxm = xm >> 32;
			unsigned lym = ym;
			unsigned hym = ym >> 32;
			unsigned long long lrm;
			unsigned long long hrm;

			lrm = DPXMULT(lxm, lym);
			hrm = DPXMULT(hxm, hym);

			{
				unsigned long long t = DPXMULT(lxm, hym);
				{
					unsigned long long at =
					    lrm + (t << 32);
					hrm += at < lrm;
					lrm = at;
				}
				hrm = hrm + (t >> 32);
			}

			{
				unsigned long long t = DPXMULT(hxm, lym);
				{
					unsigned long long at =
					    lrm + (t << 32);
					hrm += at < lrm;
					lrm = at;
				}
				hrm = hrm + (t >> 32);
			}
			rm = hrm | (lrm != 0);
		}

		/*
		 * sticky shift down to normal rounding precision
		 */
		if ((signed long long) rm < 0) {
			rm =
			    (rm >> (64 - (DP_MBITS + 1 + 3))) |
			    ((rm << (DP_MBITS + 1 + 3)) != 0);
			re++;
		} else {
			rm =
			    (rm >> (64 - (DP_MBITS + 1 + 3 + 1))) |
			    ((rm << (DP_MBITS + 1 + 3 + 1)) != 0);
		}
		assert(rm & (DP_HIDDEN_BIT << 3));
		DPNORMRET2(rs, re, rm, "mul", x, y);
	}
}

--- NEW FILE: dp_simple.c ---
/* IEEE754 floating point arithmetic
 * double precision: common utilities
 */
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.  All rights reserved.
 * http://www.algor.co.uk
 *
 * ########################################################################
 *
 *  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.,
 *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
 *
 * ########################################################################
 */


#include "ieee754dp.h"

int ieee754dp_finite(ieee754dp x)
{
	return DPBEXP(x) != DP_EMAX + 1 + DP_EBIAS;
}

ieee754dp ieee754dp_copysign(ieee754dp x, ieee754dp y)
{
	CLEARCX;
	DPSIGN(x) = DPSIGN(y);
	return x;
}


ieee754dp ieee754dp_neg(ieee754dp x)
{
	COMPXDP;

	EXPLODEXDP;
	CLEARCX;
	FLUSHXDP;

	if (xc == IEEE754_CLASS_SNAN) {
		SETCX(IEEE754_INVALID_OPERATION);
		return ieee754dp_nanxcpt(ieee754dp_indef(), "neg");
	}

	if (ieee754dp_isnan(x))	/* but not infinity */
		return ieee754dp_nanxcpt(x, "neg", x);

	/* quick fix up */
	DPSIGN(x) ^= 1;
	return x;
}


ieee754dp ieee754dp_abs(ieee754dp x)
{
	COMPXDP;

	EXPLODEXDP;
	CLEARCX;
	FLUSHXDP;

	if (xc == IEEE754_CLASS_SNAN) {
		SETCX(IEEE754_INVALID_OPERATION);
		return ieee754dp_nanxcpt(ieee754dp_indef(), "neg");
	}

	if (ieee754dp_isnan(x))	/* but not infinity */
		return ieee754dp_nanxcpt(x, "abs", x);

	/* quick fix up */
	DPSIGN(x) = 0;
	return x;
}

--- NEW FILE: dp_sqrt.c ---
/* IEEE754 floating point arithmetic
 * double precision square root
 */
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.  All rights reserved.
 * http://www.algor.co.uk
 *
 * ########################################################################
 *
 *  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.,
 *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
 *
 * ########################################################################
 */


#include "ieee754dp.h"

static const unsigned table[] = {
	0, 1204, 3062, 5746, 9193, 13348, 18162, 23592,
	29598, 36145, 43202, 50740, 58733, 67158, 75992,
	85215, 83599, 71378, 60428, 50647, 41945, 34246,
	27478, 21581, 16499, 12183, 8588, 5674, 3403,
	1742, 661, 130
};

ieee754dp ieee754dp_sqrt(ieee754dp x)
{
	struct ieee754_csr oldcsr;
	ieee754dp y, z, t;
	unsigned scalx, yh;
	COMPXDP;

	EXPLODEXDP;
	CLEARCX;
	FLUSHXDP;

	/* x == INF or NAN? */
	switch (xc) {
	case IEEE754_CLASS_QNAN:
		/* sqrt(Nan) = Nan */
		return ieee754dp_nanxcpt(x, "sqrt");
	case IEEE754_CLASS_SNAN:
		SETCX(IEEE754_INVALID_OPERATION);
		return ieee754dp_nanxcpt(ieee754dp_indef(), "sqrt");
	case IEEE754_CLASS_ZERO:
		/* sqrt(0) = 0 */
		return x;
	case IEEE754_CLASS_INF:
		if (xs) {
			/* sqrt(-Inf) = Nan */
			SETCX(IEEE754_INVALID_OPERATION);
			return ieee754dp_nanxcpt(ieee754dp_indef(), "sqrt");
		}
		/* sqrt(+Inf) = Inf */
		return x;
	case IEEE754_CLASS_DNORM:
		DPDNORMX;
		/* fall through */
	case IEEE754_CLASS_NORM:
		if (xs) {
			/* sqrt(-x) = Nan */
			SETCX(IEEE754_INVALID_OPERATION);
			return ieee754dp_nanxcpt(ieee754dp_indef(), "sqrt");
		}
		break;
	}

	/* save old csr; switch off INX enable & flag; set RN rounding */
	oldcsr = ieee754_csr;
	ieee754_csr.mx &= ~IEEE754_INEXACT;
	ieee754_csr.sx &= ~IEEE754_INEXACT;
	ieee754_csr.rm = IEEE754_RN;

	/* adjust exponent to prevent overflow */
	scalx = 0;
	if (xe > 512) {		/* x > 2**-512? */
		xe -= 512;	/* x = x / 2**512 */
		scalx += 256;
	} else if (xe < -512) {	/* x < 2**-512? */
		xe += 512;	/* x = x * 2**512 */
		scalx -= 256;
	}

	y = x = builddp(0, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);

	/* magic initial approximation to almost 8 sig. bits */
	yh = y.bits >> 32;
	yh = (yh >> 1) + 0x1ff80000;
	yh = yh - table[(yh >> 15) & 31];
	y.bits = ((unsigned long long) yh << 32) | (y.bits & 0xffffffff);

	/* Heron's rule once with correction to improve to ~18 sig. bits */
	/* t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0; */
	t = ieee754dp_div(x, y);
	y = ieee754dp_add(y, t);
	y.bits -= 0x0010000600000000LL;
	y.bits &= 0xffffffff00000000LL;

	/* triple to almost 56 sig. bits: y ~= sqrt(x) to within 1 ulp */
	/* t=y*y; z=t;  pt[n0]+=0x00100000; t+=z; z=(x-z)*y; */
	z = t = ieee754dp_mul(y, y);
	t.parts.bexp += 0x001;
	t = ieee754dp_add(t, z);
	z = ieee754dp_mul(ieee754dp_sub(x, z), y);

	/* t=z/(t+x) ;  pt[n0]+=0x00100000; y+=t; */
	t = ieee754dp_div(z, ieee754dp_add(t, x));
	t.parts.bexp += 0x001;
	y = ieee754dp_add(y, t);

	/* twiddle last bit to force y correctly rounded */

	/* set RZ, clear INEX flag */
	ieee754_csr.rm = IEEE754_RZ;
	ieee754_csr.sx &= ~IEEE754_INEXACT;

	/* t=x/y; ...chopped quotient, possibly inexact */
	t = ieee754dp_div(x, y);

	if (ieee754_csr.sx & IEEE754_INEXACT || t.bits != y.bits) {

		if (!(ieee754_csr.sx & IEEE754_INEXACT))
			/* t = t-ulp */
			t.bits -= 1;

		/* add inexact to result status */
		oldcsr.cx |= IEEE754_INEXACT;
		oldcsr.sx |= IEEE754_INEXACT;

		switch (oldcsr.rm) {
		case IEEE754_RP:
			y.bits += 1;
			/* drop through */
		case IEEE754_RN:
			t.bits += 1;
			break;
		}

		/* y=y+t; ...chopped sum */
		y = ieee754dp_add(y, t);

		/* adjust scalx for correctly rounded sqrt(x) */
		scalx -= 1;
	}

	/* py[n0]=py[n0]+scalx; ...scale back y */
	y.parts.bexp += scalx;

	/* restore rounding mode, possibly set inexact */
	ieee754_csr = oldcsr;

	return y;
}


--- NEW FILE: ieee754.c ---
/* ieee754 floating point arithmetic
 * single and double precision
 *
 * BUGS
 * not much dp done
 * doesnt generate IEEE754_INEXACT
 *
 */
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.  All rights reserved.
 * http://www.algor.co.uk
 *
 * ########################################################################
 *
 *  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.,
 *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
 *
 * ########################################################################
 */


#include "ieee754int.h"

#define DP_EBIAS	1023
#define DP_EMIN		(-1022)
#define DP_EMAX		1023

#define SP_EBIAS	127
#define SP_EMIN		(-126)
#define SP_EMAX		127

/* indexed by class */
const char *const ieee754_cname[] = {
	"Normal",
	"Zero",
	"Denormal",
	"Infinity",
	"QNaN",
	"SNaN",
};

/* the control status register 
*/
struct ieee754_csr ieee754_csr;

/* special constants
*/


#if (defined(BYTE_ORDER) && BYTE_ORDER == LITTLE_ENDIAN) || defined(__MIPSEL__)
#define SPSTR(s,b,m) {m,b,s}
#define DPSTR(s,b,mh,ml) {ml,mh,b,s}
#endif

#ifdef __MIPSEB__
#define SPSTR(s,b,m) {s,b,m}
#define DPSTR(s,b,mh,ml) {s,b,mh,ml}
#endif

const struct ieee754dp_konst __ieee754dp_spcvals[] = {
	DPSTR(0, DP_EMIN - 1 + DP_EBIAS, 0, 0),	/* + zero   */
	DPSTR(1, DP_EMIN - 1 + DP_EBIAS, 0, 0),	/* - zero   */
	DPSTR(0, DP_EBIAS, 0, 0),	/* + 1.0   */
	DPSTR(1, DP_EBIAS, 0, 0),	/* - 1.0   */
	DPSTR(0, 3 + DP_EBIAS, 0x40000, 0),	/* + 10.0   */
	DPSTR(1, 3 + DP_EBIAS, 0x40000, 0),	/* - 10.0   */
	DPSTR(0, DP_EMAX + 1 + DP_EBIAS, 0, 0),	/* + infinity */
	DPSTR(1, DP_EMAX + 1 + DP_EBIAS, 0, 0),	/* - infinity */
	DPSTR(0,DP_EMAX+1+DP_EBIAS,0x7FFFF,0xFFFFFFFF), /* + indef quiet Nan */
	DPSTR(0, DP_EMAX + DP_EBIAS, 0xFFFFF, 0xFFFFFFFF),	/* + max */
	DPSTR(1, DP_EMAX + DP_EBIAS, 0xFFFFF, 0xFFFFFFFF),	/* - max */
	DPSTR(0, DP_EMIN + DP_EBIAS, 0, 0),	/* + min normal */
	DPSTR(1, DP_EMIN + DP_EBIAS, 0, 0),	/* - min normal */
	DPSTR(0, DP_EMIN - 1 + DP_EBIAS, 0, 1),	/* + min denormal */
	DPSTR(1, DP_EMIN - 1 + DP_EBIAS, 0, 1),	/* - min denormal */
	DPSTR(0, 31 + DP_EBIAS, 0, 0),	/* + 1.0e31 */
	DPSTR(0, 63 + DP_EBIAS, 0, 0),	/* + 1.0e63 */
};

const struct ieee754sp_konst __ieee754sp_spcvals[] = {
	SPSTR(0, SP_EMIN - 1 + SP_EBIAS, 0),	/* + zero   */
	SPSTR(1, SP_EMIN - 1 + SP_EBIAS, 0),	/* - zero   */
	SPSTR(0, SP_EBIAS, 0),	/* + 1.0   */
	SPSTR(1, SP_EBIAS, 0),	/* - 1.0   */
	SPSTR(0, 3 + SP_EBIAS, 0x200000),	/* + 10.0   */
	SPSTR(1, 3 + SP_EBIAS, 0x200000),	/* - 10.0   */
	SPSTR(0, SP_EMAX + 1 + SP_EBIAS, 0),	/* + infinity */
	SPSTR(1, SP_EMAX + 1 + SP_EBIAS, 0),	/* - infinity */
	SPSTR(0,SP_EMAX+1+SP_EBIAS,0x3FFFFF),     /* + indef quiet Nan  */
	SPSTR(0, SP_EMAX + SP_EBIAS, 0x7FFFFF),	/* + max normal */
	SPSTR(1, SP_EMAX + SP_EBIAS, 0x7FFFFF),	/* - max normal */
	SPSTR(0, SP_EMIN + SP_EBIAS, 0),	/* + min normal */
	SPSTR(1, SP_EMIN + SP_EBIAS, 0),	/* - min normal */
	SPSTR(0, SP_EMIN - 1 + SP_EBIAS, 1),	/* + min denormal */
	SPSTR(1, SP_EMIN - 1 + SP_EBIAS, 1),	/* - min denormal */
	SPSTR(0, 31 + SP_EBIAS, 0),	/* + 1.0e31 */
	SPSTR(0, 63 + SP_EBIAS, 0),	/* + 1.0e63 */
};


int ieee754si_xcpt(int r, const char *op, ...)
{
	struct ieee754xctx ax;

	if (!TSTX())
		return r;
	ax.op = op;
	ax.rt = IEEE754_RT_SI;
	ax.rv.si = r;
	va_start(ax.ap, op);
	ieee754_xcpt(&ax);
	return ax.rv.si;
}

long long ieee754di_xcpt(long long r, const char *op, ...)
{
	struct ieee754xctx ax;

	if (!TSTX())
		return r;
	ax.op = op;
	ax.rt = IEEE754_RT_DI;
	ax.rv.di = r;
	va_start(ax.ap, op);
	ieee754_xcpt(&ax);
	return ax.rv.di;
}

--- NEW FILE: ieee754.h ---
/* single and double precision fp ops
 * missing extended precision.
*/
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.  All rights reserved.
 * http://www.algor.co.uk
 *
 * ########################################################################
 *
 *  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.,
 *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
 *
 * ########################################################################
 */

/**************************************************************************
 *  Nov 7, 2000
 *  Modification to allow integration with Linux kernel 
 *
 *  Kevin D. Kissell, ke...@mi... and Carsten Langgard, car...@mi...
 *  Copyright (C) 2000 MIPS Technologies, Inc. All rights reserved.
 *************************************************************************/

#ifdef __KERNEL__
/* Going from Algorithmics to Linux native environment, add this */
#include <linux/types.h>

/* 
 * Not very pretty, but the Linux kernel's normal va_list definition 
 * does not allow it to be used as a structure element, as it is here.
 */
#ifndef _STDARG_H
#include <stdarg.h>
#endif

#else

/* Note that __KERNEL__ is taken to mean Linux kernel */

#if #system(OpenBSD)
#include <machine/types.h>
#endif
#include <machine/endian.h>

#endif				/* __KERNEL__ */

#if (defined(BYTE_ORDER) && BYTE_ORDER == LITTLE_ENDIAN) || defined(__MIPSEL__)
struct ieee754dp_konst {
	unsigned mantlo:32;
	unsigned manthi:20;
	unsigned bexp:11;
	unsigned sign:1;
};
struct ieee754sp_konst {
	unsigned mant:23;
	unsigned bexp:8;
	unsigned sign:1;
};

typedef union _ieee754dp {
	struct ieee754dp_konst oparts;
	struct {
		unsigned long long mant:52;
		unsigned int bexp:11;
		unsigned int sign:1;
	} parts;
	unsigned long long bits;
	double d;
} ieee754dp;

typedef union _ieee754sp {
	struct ieee754sp_konst parts;
	float f;
	unsigned long bits;
} ieee754sp;
#endif

#if (defined(BYTE_ORDER) && BYTE_ORDER == BIG_ENDIAN) || defined(__MIPSEB__)
struct ieee754dp_konst {
	unsigned sign:1;
	unsigned bexp:11;
	unsigned manthi:20;
	unsigned mantlo:32;
};
typedef union _ieee754dp {
	struct ieee754dp_konst oparts;
	struct {
		unsigned int sign:1;
		unsigned int bexp:11;
		unsigned long long mant:52;
	} parts;
	double d;
	unsigned long long bits;
} ieee754dp;

struct ieee754sp_konst {
	unsigned sign:1;
	unsigned bexp:8;
	unsigned mant:23;
};

typedef union _ieee754sp {
	struct ieee754sp_konst parts;
	float f;
	unsigned long bits;
} ieee754sp;
#endif

/*
 * single precision (often aka float)
*/
int ieee754sp_finite(ieee754sp x);
int ieee754sp_class(ieee754sp x);

ieee754sp ieee754sp_abs(ieee754sp x);
ieee754sp ieee754sp_neg(ieee754sp x);
ieee754sp ieee754sp_scalb(ieee754sp x, int);
ieee754sp ieee754sp_logb(ieee754sp x);

/* x with sign of y */
ieee754sp ieee754sp_copysign(ieee754sp x, ieee754sp y);

ieee754sp ieee754sp_add(ieee754sp x, ieee754sp y);
ieee754sp ieee754sp_sub(ieee754sp x, ieee754sp y);
ieee754sp ieee754sp_mul(ieee754sp x, ieee754sp y);
ieee754sp ieee754sp_div(ieee754sp x, ieee754sp y);

ieee754sp ieee754sp_fint(int x);
ieee754sp ieee754sp_funs(unsigned x);
ieee754sp ieee754sp_flong(long long x);
ieee754sp ieee754sp_fulong(unsigned long long x);
ieee754sp ieee754sp_fdp(ieee754dp x);

int ieee754sp_tint(ieee754sp x);
unsigned int ieee754sp_tuns(ieee754sp x);
long long ieee754sp_tlong(ieee754sp x);
unsigned long long ieee754sp_tulong(ieee754sp x);

int ieee754sp_cmp(ieee754sp x, ieee754sp y, int cop);
/*
 * basic sp math
 */
ieee754sp ieee754sp_modf(ieee754sp x, ieee754sp * ip);
ieee754sp ieee754sp_frexp(ieee754sp x, int *exp);
ieee754sp ieee754sp_ldexp(ieee754sp x, int exp);

ieee754sp ieee754sp_ceil(ieee754sp x);
ieee754sp ieee754sp_floor(ieee754sp x);
ieee754sp ieee754sp_trunc(ieee754sp x);

ieee754sp ieee754sp_sqrt(ieee754sp x);

/*
 * double precision (often aka double)
*/
int ieee754dp_finite(ieee754dp x);
int ieee754dp_class(ieee754dp x);

/* x with sign of y */
ieee754dp ieee754dp_copysign(ieee754dp x, ieee754dp y);

ieee754dp ieee754dp_add(ieee754dp x, ieee754dp y);
ieee754dp ieee754dp_sub(ieee754dp x, ieee754dp y);
ieee754dp ieee754dp_mul(ieee754dp x, ieee754dp y);
ieee754dp ieee754dp_div(ieee754dp x, ieee754dp y);

ieee754dp ieee754dp_abs(ieee754dp x);
ieee754dp ieee754dp_neg(ieee754dp x);
ieee754dp ieee754dp_scalb(ieee754dp x, int);

/* return exponent as integer in floating point format 
 */
ieee754dp ieee754dp_logb(ieee754dp x);

ieee754dp ieee754dp_fint(int x);
ieee754dp ieee754dp_funs(unsigned x);
ieee754dp ieee754dp_flong(long long x);
ieee754dp ieee754dp_fulong(unsigned long long x);
ieee754dp ieee754dp_fsp(ieee754sp x);

ieee754dp ieee754dp_ceil(ieee754dp x);
ieee754dp ieee754dp_floor(ieee754dp x);
ieee754dp ieee754dp_trunc(ieee754dp x);

int ieee754dp_tint(ieee754dp x);
unsigned int ieee754dp_tuns(ieee754dp x);
long long ieee754dp_tlong(ieee754dp x);
unsigned long long ieee754dp_tulong(ieee754dp x);

int ieee754dp_cmp(ieee754dp x, ieee754dp y, int cop);
/*
 * basic sp math
 */
ieee754dp ieee754dp_modf(ieee754dp x, ieee754dp * ip);
ieee754dp ieee754dp_frexp(ieee754dp x, int *exp);
ieee754dp ieee754dp_ldexp(ieee754dp x, int exp);

ieee754dp ieee754dp_ceil(ieee754dp x);
ieee754dp ieee754dp_floor(ieee754dp x);
ieee754dp ieee754dp_trunc(ieee754dp x);

ieee754dp ieee754dp_sqrt(ieee754dp x);



/* 5 types of floating point number 
*/
#define IEEE754_CLASS_NORM	0x00
#define IEEE754_CLASS_ZERO	0x01
#define IEEE754_CLASS_DNORM	0x02
#define IEEE754_CLASS_INF	0x03
#define IEEE754_CLASS_SNAN	0x04
#define IEEE754_CLASS_QNAN	0x05
extern const char *const ieee754_cname[];

/* exception numbers */
#define IEEE754_INEXACT			0x01
#define IEEE754_UNDERFLOW		0x02
#define IEEE754_OVERFLOW		0x04
#define IEEE754_ZERO_DIVIDE		0x08
#define IEEE754_INVALID_OPERATION	0x10

/* cmp operators
*/
#define IEEE754_CLT	0x01
#define IEEE754_CEQ	0x02
#define IEEE754_CGT	0x04
#define IEEE754_CUN	0x08

/* rounding mode 
*/
#define IEEE754_RN	0	/* round to nearest */
#define IEEE754_RZ	1	/* round toward zero  */
#define IEEE754_RD	2	/* round toward -Infinity */
#define IEEE754_RU	3	/* round toward +Infinity */

/* other naming */
#define IEEE754_RM	IEEE754_RD
#define IEEE754_RP	IEEE754_RU

/* "normal" comparisons
*/
static __inline int ieee754sp_eq(ieee754sp x, ieee754sp y)
{
	return ieee754sp_cmp(x, y, IEEE754_CEQ);
}

static __inline int ieee754sp_ne(ieee754sp x, ieee754sp y)
{
	return ieee754sp_cmp(x, y,
			     IEEE754_CLT | IEEE754_CGT | IEEE754_CUN);
}

static __inline int ieee754sp_lt(ieee754sp x, ieee754sp y)
{
	return ieee754sp_cmp(x, y, IEEE754_CLT);
}

static __inline int ieee754sp_le(ieee754sp x, ieee754sp y)
{
	return ieee754sp_cmp(x, y, IEEE754_CLT | IEEE754_CEQ);
}

static __inline int ieee754sp_gt(ieee754sp x, ieee754sp y)
{
	return ieee754sp_cmp(x, y, IEEE754_CGT);
}


static __inline int ieee754sp_ge(ieee754sp x, ieee754sp y)
{
	return ieee754sp_cmp(x, y, IEEE754_CGT | IEEE754_CEQ);
}

static __inline int ieee754dp_eq(ieee754dp x, ieee754dp y)
{
	return ieee754dp_cmp(x, y, IEEE754_CEQ);
}

static __inline int ieee754dp_ne(ieee754dp x, ieee754dp y)
{
	return ieee754dp_cmp(x, y,
			     IEEE754_CLT | IEEE754_CGT | IEEE754_CUN);
}

static __inline int ieee754dp_lt(ieee754dp x, ieee754dp y)
{
	return ieee754dp_cmp(x, y, IEEE754_CLT);
}

static __inline int ieee754dp_le(ieee754dp x, ieee754dp y)
{
	return ieee754dp_cmp(x, y, IEEE754_CLT | IEEE754_CEQ);
}

static __inline int ieee754dp_gt(ieee754dp x, ieee754dp y)
{
	return ieee754dp_cmp(x, y, IEEE754_CGT);
}

static __inline int ieee754dp_ge(ieee754dp x, ieee754dp y)
{
	return ieee754dp_cmp(x, y, IEEE754_CGT | IEEE754_CEQ);
}


/* like strtod
*/
ieee754dp ieee754dp_fstr(const char *s, char **endp);
char *ieee754dp_tstr(ieee754dp x, int prec, int fmt, int af);


/* the control status register 
*/
struct ieee754_csr {
	unsigned pad:13;
	unsigned nod:1;		/* set 1 for no denormalised numbers */
	unsigned cx:5;		/* exceptions this operation */
	unsigned mx:5;		/* exception enable  mask */
	unsigned sx:5;		/* exceptions total */
	unsigned rm:2;		/* current rounding mode */
};
extern struct ieee754_csr ieee754_csr;

static __inline unsigned ieee754_getrm(void)
{
	return (ieee754_csr.rm);
}
static __inline unsigned ieee754_setrm(unsigned rm)
{
	return (ieee754_csr.rm = rm);
}

/*
 * get current exceptions
 */
static __inline unsigned ieee754_getcx(void)
{
	return (ieee754_csr.cx);
}

/* test for current exception condition 
 */
static __inline int ieee754_cxtest(unsigned n)
{
	return (ieee754_csr.cx & n);
}

/*
 * get sticky exceptions
 */
static __inline unsigned ieee754_getsx(void)
{
	return (ieee754_csr.sx);
}

/* clear sticky conditions
*/
static __inline unsigned ieee754_clrsx(void)
{
	return (ieee754_csr.sx = 0);
}

/* test for sticky exception condition 
 */
static __inline int ieee754_sxtest(unsigned n)
{
	return (ieee754_csr.sx & n);
}

/* debugging */
ieee754sp ieee754sp_dump(char *s, ieee754sp x);
ieee754dp ieee754dp_dump(char *s, ieee754dp x);

#define IEEE754_SPCVAL_PZERO	0
#define IEEE754_SPCVAL_NZERO	1
#define IEEE754_SPCVAL_PONE	2
#define IEEE754_SPCVAL_NONE	3
#define IEEE754_SPCVAL_PTEN	4
#define IEEE754_SPCVAL_NTEN	5
#define IEEE754_SPCVAL_PINFINITY	6
#define IEEE754_SPCVAL_NINFINITY	7
#define IEEE754_SPCVAL_INDEF	8
#define IEEE754_SPCVAL_PMAX	9	/* +max norm */
#define IEEE754_SPCVAL_NMAX	10	/* -max norm */
#define IEEE754_SPCVAL_PMIN	11	/* +min norm */
#define IEEE754_SPCVAL_NMIN	12	/* +min norm */
#define IEEE754_SPCVAL_PMIND	13	/* +min denorm */
#define IEEE754_SPCVAL_NMIND	14	/* +min denorm */
#define IEEE754_SPCVAL_P1E31	15	/* + 1.0e31 */
#define IEEE754_SPCVAL_P1E63	16	/* + 1.0e63 */

extern const struct ieee754dp_konst __ieee754dp_spcvals[];
extern const struct ieee754sp_konst __ieee754sp_spcvals[];
#define ieee754dp_spcvals ((const ieee754dp *)__ieee754dp_spcvals)
#define ieee754sp_spcvals ((const ieee754sp *)__ieee754sp_spcvals)

/* return infinity with given sign
*/
#define ieee754dp_inf(sn)	\
  (ieee754dp_spcvals[IEEE754_SPCVAL_PINFINITY+(sn)])
#define ieee754dp_zero(sn) \
  (ieee754dp_spcvals[IEEE754_SPCVAL_PZERO+(sn)])
#define ieee754dp_one(sn) \
  (ieee754dp_spcvals[IEEE754_SPCVAL_PONE+(sn)])
#define ieee754dp_ten(sn) \
  (ieee754dp_spcvals[IEEE754_SPCVAL_PTEN+(sn)])
#define ieee754dp_indef() \
  (ieee754dp_spcvals[IEEE754_SPCVAL_INDEF])
#define ieee754dp_max(sn) \
  (ieee754dp_spcvals[IEEE754_SPCVAL_PMAX+(sn)])
#define ieee754dp_min(sn) \
  (ieee754dp_spcvals[IEEE754_SPCVAL_PMIN+(sn)])
#define ieee754dp_mind(sn) \
  (ieee754dp_spcvals[IEEE754_SPCVAL_PMIND+(sn)])
#define ieee754dp_1e31() \
  (ieee754dp_spcvals[IEEE754_SPCVAL_P1E31])
#define ieee754dp_1e63() \
  (ieee754dp_spcvals[IEEE754_SPCVAL_P1E63])

#define ieee754sp_inf(sn) \
  (ieee754sp_spcvals[IEEE754_SPCVAL_PINFINITY+(sn)])
#define ieee754sp_zero(sn) \
  (ieee754sp_spcvals[IEEE754_SPCVAL_PZERO+(sn)])
#define ieee754sp_one(sn) \
  (ieee754sp_spcvals[IEEE754_SPCVAL_PONE+(sn)])
#define ieee754sp_ten(sn) \
  (ieee754sp_spcvals[IEEE754_SPCVAL_PTEN+(sn)])
#define ieee754sp_indef() \
  (ieee754sp_spcvals[IEEE754_SPCVAL_INDEF])
#define ieee754sp_max(sn) \
  (ieee754sp_spcvals[IEEE754_SPCVAL_PMAX+(sn)])
#define ieee754sp_min(sn) \
  (ieee754sp_spcvals[IEEE754_SPCVAL_PMIN+(sn)])
#define ieee754sp_mind(sn) \
  (ieee754sp_spcvals[IEEE754_SPCVAL_PMIND+(sn)])
#define ieee754sp_1e31() \
  (ieee754sp_spcvals[IEEE754_SPCVAL_P1E31])
#define ieee754sp_1e63() \
  (ieee754sp_spcvals[IEEE754_SPCVAL_P1E63])

/* indefinite integer value 
*/
#define ieee754si_indef()	INT_MAX
#ifdef LONG_LONG_MAX
#define ieee754di_indef()	LONG_LONG_MAX
#else
#define ieee754di_indef()	((long long)(~0ULL>>1))
#endif

/* IEEE exception context, passed to handler */
struct ieee754xctx {
	const char *op;		/* operation name */
	int rt;			/* result type */
	union {
		ieee754sp sp;	/* single precision */
		ieee754dp dp;	/* double precision */
#ifdef IEEE854_XP
		ieee754xp xp;	/* extended precision */
#endif
		int si;		/* standard signed integer (32bits) */
		long long di;	/* extended signed integer (64bits) */
	} rv;			/* default result format implied by op */
	va_list ap;
};

/* result types for xctx.rt */
#define IEEE754_RT_SP	0
#define IEEE754_RT_DP	1
#define IEEE754_RT_XP	2
#define IEEE754_RT_SI	3
#define IEEE754_RT_DI	4

extern void ieee754_xcpt(struct ieee754xctx *xcp);

/* compat */
#define ieee754dp_fix(x)	ieee754dp_tint(x)
#define ieee754sp_fix(x)	ieee754sp_tint(x)

--- NEW FILE: ieee754dp.c ---
/* IEEE754 floating point arithmetic
 * double precision: common utilities
 */
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.  All rights reserved.
 * http://www.algor.co.uk
 *
 * ########################################################################
 *
 *  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.,
 *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
 *
 * ########################################################################
 */


#include "ieee754dp.h"

int ieee754dp_class(ieee754dp x)
{
	COMPXDP;
	EXPLODEXDP;
	return xc;
}

int ieee754dp_isnan(ieee754dp x)
{
	return ieee754dp_class(x) >= IEEE754_CLASS_SNAN;
}

int ieee754dp_issnan(ieee754dp x)
{
	assert(ieee754dp_isnan(x));
	return ((DPMANT(x) & DP_MBIT(DP_MBITS-1)) == DP_MBIT(DP_MBITS-1));
}


ieee754dp ieee754dp_xcpt(ieee754dp r, const char *op, ...)
{
	struct ieee754xctx ax;
	if (!TSTX())
		return r;

	ax.op = op;
	ax.rt = IEEE754_RT_DP;
	ax.rv.dp = r;
	va_start(ax.ap, op);
	ieee754_xcpt(&ax);
	return ax.rv.dp;
}

ieee754dp ieee754dp_nanxcpt(ieee754dp r, const char *op, ...)
{
	struct ieee754xctx ax;

	assert(ieee754dp_isnan(r));

	if (!ieee754dp_issnan(r))	/* QNAN does not cause invalid op !! */
		return r;

	if (!SETCX(IEEE754_INVALID_OPERATION)) {
		/* not enabled convert to a quiet NaN */
		DPMANT(r) &= (~DP_MBIT(DP_MBITS-1));
		if (ieee754dp_isnan(r))
			return r;
		else
			return ieee754dp_indef();
	}

	ax.op = op;
	ax.rt = 0;
	ax.rv.dp = r;
	va_start(ax.ap, op);
	ieee754_xcpt(&ax);
	return ax.rv.dp;
}

ieee754dp ieee754dp_bestnan(ieee754dp x, ieee754dp y)
{
	assert(ieee754dp_isnan(x));
	assert(ieee754dp_isnan(y));

	if (DPMANT(x) > DPMANT(y))
		return x;
	else
		return y;
}


/* generate a normal/denormal number with over,under handeling
 * sn is sign
 * xe is an unbiased exponent
 * xm is 3bit extended precision value.
 */
ieee754dp ieee754dp_format(int sn, int xe, unsigned long long xm)
{
	assert(xm);		/* we dont gen exact zeros (probably should) */

	assert((xm >> (DP_MBITS + 1 + 3)) == 0);	/* no execess */
	assert(xm & (DP_HIDDEN_BIT << 3));

	if (xe < DP_EMIN) {
		/* strip lower bits */
		int es = DP_EMIN - xe;

		if (ieee754_csr.nod) {
			SETCX(IEEE754_UNDERFLOW);
			SETCX(IEEE754_INEXACT);

			switch(ieee754_csr.rm) {
			case IEEE754_RN:
				return ieee754dp_zero(sn);
			case IEEE754_RZ:
				return ieee754dp_zero(sn);
			case IEEE754_RU:    /* toward +Infinity */
				if(sn == 0)
					return ieee754dp_min(0);
				else
					return ieee754dp_zero(1);
			case IEEE754_RD:    /* toward -Infinity */
				if(sn == 0)
					return ieee754dp_zero(0);
				else
					return ieee754dp_min(1);
			}
		}

		/* sticky right shift es bits 
		 */
		xm = XDPSRS(xm, es);
		xe += es;

		assert((xm & (DP_HIDDEN_BIT << 3)) == 0);
		assert(xe == DP_EMIN);
	}
	if (xm & (DP_MBIT(3) - 1)) {
		SETCX(IEEE754_INEXACT);
		/* inexact must round of 3 bits 
		 */
		switch (ieee754_csr.rm) {
		case IEEE754_RZ:
			break;
		case IEEE754_RN:
			xm += 0x3 + ((xm >> 3) & 1);
			/* xm += (xm&0x8)?0x4:0x3 */
			break;
		case IEEE754_RU:	/* toward +Infinity */
			if (!sn)	/* ?? */
				xm += 0x8;
			break;
		case IEEE754_RD:	/* toward -Infinity */
			if (sn)	/* ?? */
				xm += 0x8;
			break;
		}
		/* adjust exponent for rounding add overflowing 
		 */
		if (xm >> (DP_MBITS + 3 + 1)) {	/* add causes mantissa overflow */
			xm >>= 1;
			xe++;
		}
	}
	/* strip grs bits */
	xm >>= 3;

	assert((xm >> (DP_MBITS + 1)) == 0);	/* no execess */
	assert(xe >= DP_EMIN);

	if (xe > DP_EMAX) {
		SETCX(IEEE754_OVERFLOW);
		SETCX(IEEE754_INEXACT);
		/* -O can be table indexed by (rm,sn) */
		switch (ieee754_csr.rm) {
		case IEEE754_RN:
			return ieee754dp_inf(sn);
		case IEEE754_RZ:
			return ieee754dp_max(sn);
		case IEEE754_RU:	/* toward +Infinity */
			if (sn == 0)
				return ieee754dp_inf(0);
			else
				return ieee754dp_max(1);
		case IEEE754_RD:	/* toward -Infinity */
			if (sn == 0)
				return ieee754dp_max(0);
			else
				return ieee754dp_inf(1);
		}
	}
	/* gen norm/denorm/zero */

	if ((xm & DP_HIDDEN_BIT) == 0) {
		/* we underflow (tiny/zero) */
		assert(xe == DP_EMIN);
		SETCX(IEEE754_UNDERFLOW);
		return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm);
	} else {
		assert((xm >> (DP_MBITS + 1)) == 0);	/* no execess */
		assert(xm & DP_HIDDEN_BIT);

		return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
	}
}

--- NEW FILE: ieee754int.h ---
/* 
 * IEEE754 floating point
 * common internal header file
 */
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.  All rights reserved.
 * http://www.algor.co.uk
 *
 * ########################################################################
 *
 *  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.,
 *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
 *
 * ########################################################################
 */


#include "ieee754.h"

#define DP_EBIAS	1023
#define DP_EMIN		(-1022)
#define DP_EMAX		1023
#define DP_MBITS	52

#define SP_EBIAS	127
#define SP_EMIN		(-126)
#define SP_EMAX		127
#define SP_MBITS	23

#define DP_MBIT(x)	((unsigned long long)1 << (x))
#define DP_HIDDEN_BIT	DP_MBIT(DP_MBITS)
#define DP_SIGN_BIT	DP_MBIT(63)

#define SP_MBIT(x)	((unsigned long)1 << (x))
#define SP_HIDDEN_BIT	SP_MBIT(SP_MBITS)
#define SP_SIGN_BIT	SP_MBIT(31)


#define SPSIGN(sp)	(sp.parts.sign)
#define SPBEXP(sp)	(sp.parts.bexp)
#define SPMANT(sp)	(sp.parts.mant)

#define DPSIGN(dp)	(dp.parts.sign)
#define DPBEXP(dp)	(dp.parts.bexp)
#define DPMANT(dp)	(dp.parts.mant)

#define CLPAIR(x,y)	((x)*6+(y))

#define CLEARCX	\
  (ieee754_csr.cx = 0)

#define SETCX(x) \
  (ieee754_csr.cx |= (x),ieee754_csr.sx |= (x),ieee754_csr.mx & (x))

#define TSTX()	\
	(ieee754_csr.cx & ieee754_csr.mx)


#define COMPXSP \
  unsigned xm; int xe; int xs; int xc

#define COMPYSP \
  unsigned ym; int ye; int ys; int yc

#define EXPLODESP(v,vc,vs,ve,vm) \
{\
    vs = SPSIGN(v);\
    ve = SPBEXP(v);\
    vm = SPMANT(v);\
    if(ve == SP_EMAX+1+SP_EBIAS){\
	if(vm == 0)\
	  vc = IEEE754_CLASS_INF;\
	else if(vm & SP_MBIT(SP_MBITS-1)) \
	  vc = IEEE754_CLASS_SNAN;\
	else \
	  vc = IEEE754_CLASS_QNAN;\
    } else if(ve == SP_EMIN-1+SP_EBIAS) {\
	if(vm) {\
	    ve = SP_EMIN;\
	    vc = IEEE754_CLASS_DNORM;\
	} else\
	  vc = IEEE754_CLASS_ZERO;\
    } else {\
	ve -= SP_EBIAS;\
	vm |= SP_HIDDEN_BIT;\
	vc = IEEE754_CLASS_NORM;\
    }\
}
#define EXPLODEXSP EXPLODESP(x,xc,xs,xe,xm)
#define EXPLODEYSP EXPLODESP(y,yc,ys,ye,ym)


#define COMPXDP \
unsigned long long xm; int xe; int xs; int xc

#define COMPYDP \
unsigned long long ym; int ye; int ys; int yc

#define EXPLODEDP(v,vc,vs,ve,vm) \
{\
    vm = DPMANT(v);\
    vs = DPSIGN(v);\
    ve = DPBEXP(v);\
    if(ve == DP_EMAX+1+DP_EBIAS){\
	if(vm == 0)\
	  vc = IEEE754_CLASS_INF;\
	else if(vm & DP_MBIT(DP_MBITS-1)) \
	  vc = IEEE754_CLASS_SNAN;\
	else \
	  vc = IEEE754_CLASS_QNAN;\
    } else if(ve == DP_EMIN-1+DP_EBIAS) {\
	if(vm) {\
	    ve = DP_EMIN;\
	    vc = IEEE754_CLASS_DNORM;\
	} else\
	  vc = IEEE754_CLASS_ZERO;\
    } else {\
	ve -= DP_EBIAS;\
	vm |= DP_HIDDEN_BIT;\
	vc = IEEE754_CLASS_NORM;\
    }\
}
#define EXPLODEXDP EXPLODEDP(x,xc,xs,xe,xm)
#define EXPLODEYDP EXPLODEDP(y,yc,ys,ye,ym)

#define FLUSHDP(v,vc,vs,ve,vm) \
	if(vc==IEEE754_CLASS_DNORM) {\
	    if(ieee754_csr.nod) {\
		SETCX(IEEE754_INEXACT);\
		vc = IEEE754_CLASS_ZERO;\
		ve = DP_EMIN-1+DP_EBIAS;\
		vm = 0;\
		v = ieee754dp_zero(vs);\
	    }\
	}

#define FLUSHSP(v,vc,vs,ve,vm) \
	if(vc==IEEE754_CLASS_DNORM) {\
	    if(ieee754_csr.nod) {\
		SETCX(IEEE754_INEXACT);\
		vc = IEEE754_CLASS_ZERO;\
		ve = SP_EMIN-1+SP_EBIAS;\
		vm = 0;\
		v = ieee754sp_zero(vs);\
	    }\
	}

#define FLUSHXDP FLUSHDP(x,xc,xs,xe,xm)
#define FLUSHYDP FLUSHDP(y,yc,ys,ye,ym)
#define FLUSHXSP FLUSHSP(x,xc,xs,xe,xm)
#define FLUSHYSP FLUSHSP(y,yc,ys,ye,ym)

--- NEW FILE: ieee754sp.c ---
/* IEEE754 floating point arithmetic
 * single precision
 */
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.  All rights reserved.
 * http://www.algor.co.uk
 *
 * ########################################################################
 *
 *  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.,
 *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
 *
 * ########################################################################
 */


#include "ieee754sp.h"

int ieee754sp_class(ieee754sp x)
{
	COMPXSP;
	EXPLODEXSP;
	return xc;
}

int ieee754sp_isnan(ieee754sp x)
{
	return ieee754sp_class(x) >= IEEE754_CLASS_SNAN;
}

int ieee754sp_issnan(ieee754sp x)
{
	assert(ieee754sp_isnan(x));
	return (SPMANT(x) & SP_MBIT(SP_MBITS-1));
}


ieee754sp ieee754sp_xcpt(ieee754sp r, const char *op, ...)
{
	struct ieee754xctx ax;

	if (!TSTX())
		return r;

	ax.op = op;
	ax.rt = IEEE754_RT_SP;
	ax.rv.sp = r;
	va_start(ax.ap, op);
	ieee754_xcpt(&ax);
	return ax.rv.sp;
}

ieee754sp ieee754sp_nanxcpt(ieee754sp r, const char *op, ...)
{
	struct ieee754xctx ax;

	assert(ieee754sp_isnan(r));

	if (!ieee754sp_issnan(r))	/* QNAN does not cause invalid op !! */
		return r;

	if (!SETCX(IEEE754_INVALID_OPERATION)) {
		/* not enabled convert to a quiet NaN */
		SPMANT(r) &= (~SP_MBIT(SP_MBITS-1));
		if (ieee754sp_isnan(r))
			return r;
		else
			return ieee754sp_indef();
	}

	ax.op = op;
	ax.rt = 0;
	ax.rv.sp = r;
	va_start(ax.ap, op);
	ieee754_xcpt(&ax);
	return ax.rv.sp;
}

ieee754sp ieee754sp_bestnan(ieee754sp x, ieee754sp y)
{
	assert(ieee754sp_isnan(x));
	assert(ieee754sp_isnan(y));

	if (SPMANT(x) > SPMANT(y))
		return x;
	else
		return y;
}


/* generate a normal/denormal number with over,under handeling
 * sn is sign
 * xe is an unbiased exponent
 * xm is 3bit extended precision value.
 */
ieee754sp ieee754sp_format(int sn, int xe, unsigned xm)
{
	assert(xm);		/* we dont gen exact zeros (probably should) */

	assert((xm >> (SP_MBITS + 1 + 3)) == 0);	/* no execess */
	assert(xm & (SP_HIDDEN_BIT << 3));

	if (xe < SP_EMIN) {
		/* strip lower bits */
		int es = SP_EMIN - xe;

		if (ieee754_csr.nod) {
		    SETCX(IEEE754_UNDERFLOW);
		    SETCX(IEEE754_INEXACT);

		    switch(ieee754_csr.rm) {
			case IEEE754_RN:
			    return ieee754sp_zero(sn);
			case IEEE754_RZ:
			    return ieee754sp_zero(sn);
			case IEEE754_RU:      /* toward +Infinity */
			    if(sn == 0)
				return ieee754sp_min(0);
			    else
				return ieee754sp_zero(1);
			case IEEE754_RD:      /* toward -Infinity */
			    if(sn == 0)
				return ieee754sp_zero(0);
			    else
				return ieee754sp_min(1);
		    }
		}

		/* sticky right shift es bits 
		 */
		SPXSRSXn(es);

		assert((xm & (SP_HIDDEN_BIT << 3)) == 0);
		assert(xe == SP_EMIN);
	}
	if (xm & (SP_MBIT(3) - 1)) {
		SETCX(IEEE754_INEXACT);
		/* inexact must round of 3 bits 
		 */
		switch (ieee754_csr.rm) {
		case IEEE754_RZ:
			break;
		case IEEE754_RN:
			xm += 0x3 + ((xm >> 3) & 1);
			/* xm += (xm&0x8)?0x4:0x3 */
			break;
		case IEEE754_RU:	/* toward +Infinity */
			if (!sn)	/* ?? */
				xm += 0x8;
			break;
		case IEEE754_RD:	/* toward -Infinity */
			if (sn)	/* ?? */
				xm += 0x8;
			break;
		}
		/* adjust exponent for rounding add overflowing 
		 */
		if (xm >> (SP_MBITS + 1 + 3)) {	/* add causes mantissa overflow */
			xm >>= 1;
			xe++;
		}
	}
	/* strip grs bits */
	xm >>= 3;

	assert((xm >> (SP_MBITS + 1)) == 0);	/* no execess */
	assert(xe >= SP_EMIN);

	if (xe > SP_EMAX) {
		SETCX(IEEE754_OVERFLOW);
		SETCX(IEEE754_INEXACT);
		/* -O can be table indexed by (rm,sn) */
		switch (ieee754_csr.rm) {
		case IEEE754_RN:
			return ieee754sp_inf(sn);
		case IEEE754_RZ:
			return ieee754sp_max(sn);
		case IEEE754_RU:	/* toward +Infinity */
			if (sn == 0)
				return ieee754sp_inf(0);
			else
				return ieee754sp_max(1);
		case IEEE754_RD:	/* toward -Infinity */
			if (sn == 0)
				return ieee754sp_max(0);
			else
				return ieee754sp_inf(1);
		}
	}
	/* gen norm/denorm/zero */

	if ((xm & SP_HIDDEN_BIT) == 0) {
		/* we underflow (tiny/zero) */
		assert(xe == SP_EMIN);
		SETCX(IEEE754_UNDERFLOW);
		return buildsp(sn, SP_EMIN - 1 + SP_EBIAS, xm);
	} else {
		assert((xm >> (SP_MBITS + 1)) == 0);	/* no execess */
		assert(xm & SP_HIDDEN_BIT);

		return buildsp(sn, xe + SP_EBIAS, xm & ~SP_HIDDEN_BIT);
	}
}

--- NEW FILE: kernel_linkage.c ---
/**************************************************************************
 *
 *  arch/mips/math_emu/kernel_linkage.c
 *
 *  Kevin D. Kissell, kevink@mips and Carsten Langgaard, car...@mi...
 *  Copyright (C) 2000 MIPS Technologies, Inc.  All rights reserved.
 *
 * ########################################################################
 *
 *  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.,
 *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
 *
 *************************************************************************/
/*
 * Routines corresponding to Linux kernel FP context
 * manipulation primitives for the Algorithmics MIPS
 * FPU Emulator
 */

#include <linux/sched.h>
#include <asm/processor.h>
#include <asm/signal.h>
#include <asm/uaccess.h>

#include <asm/fpu_emulator.h>

extern struct mips_fpu_emulator_private fpuemuprivate;

#define SIGNALLING_NAN 0x7ff800007ff80000LL

void fpu_emulator_init_fpu(void)
{
	static int first = 1;
	int i;
 
	if (first) {
		first = 0;
		printk("Algorithmics/MIPS FPU Emulator v1.5\n");
	}

	current->thread.fpu.soft.sr = 0;
	for (i = 0; i < 32; i++) {
		current->thread.fpu.soft.regs[i] = SIGNALLING_NAN;
	}
}


/*
 * Emulator context save/restore to/from a signal context
 * presumed to be on the user stack, and therefore accessed
 * with appropriate macros from uaccess.h
 */

int fpu_emulator_save_context(struct sigcontext *sc)
{
	int i;
	int err = 0;

	for (i = 0; i < 32; i++) {
		err |=
		    __put_user(current->thread.fpu.soft.regs[i],
			       &sc->sc_fpregs[i]);
	}
	err |= __put_user(current->thread.fpu.soft.sr, &sc->sc_fpc_csr);
	err |= __put_user(fpuemuprivate.eir, &sc->sc_fpc_eir);

	return err;
}

int fpu_emulator_restore_context(struct sigcontext *sc)
{
	int i;
	int err = 0;

	for (i = 0; i < 32; i++) {
		err |=
		    __get_user(current->thread.fpu.soft.regs[i],
			       &sc->sc_fpregs[i]);
	}
	err |= __get_user(current->thread.fpu.soft.sr, &sc->sc_fpc_csr);
	err |= __get_user(fpuemuprivate.eir, &sc->sc_fpc_eir);

	return err;
}


--- NEW FILE: sp_add.c ---
/* IEEE754 floating point arithmetic
 * single precision
 */
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.  All rights reserved.
 * http://www.algor.co.uk
 *
 * ########################################################################
 *
 *  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.,
 *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
 *
 * ########################################################################
 */


#include "ieee754sp.h"

ieee754sp ieee754sp_add(ieee754sp x, ieee754sp y)
{
	COMPXSP;
	COMPYSP;

	EXPLODEXSP;
	EXPLODEYSP;

	CLEARCX;

	FLUSHXSP;
	FLUSHYSP;

	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):
		SETCX(IEEE754_INVALID_OPERATION);
		return ieee754sp_nanxcpt(ieee754sp_indef(), "add", x, y);

	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):
	c...
 
[truncated message content]