lk/linux-kernel-test
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Linux-2.6.12-rc2

Browse Source 
Initial git repository build. I'm not bothering with the full history,
even though we have it. We can create a separate "historical" git
archive of that later if we want to, and in the meantime it's about
3.2GB when imported into git - space that would just make the early
git days unnecessarily complicated, when we don't have a lot of good
infrastructure for it.

Let it rip!
This commit is contained in:
Linus Torvalds
2005-04-16 15:20:36 -07:00
commit 1da177e4c3
17291 changed files with 6718755 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
arch/mips/math-emu/Makefile Normal file
View File

@@ -0,0 +1,11 @@
#
# Makefile for the Linux/MIPS kernel FPU emulation.
#
obj-y := cp1emu.o ieee754m.o ieee754d.o ieee754dp.o ieee754sp.o ieee754.o \
ieee754xcpt.o dp_frexp.o dp_modf.o dp_div.o dp_mul.o dp_sub.o \
dp_add.o dp_fsp.o dp_cmp.o dp_logb.o dp_scalb.o dp_simple.o \
dp_tint.o dp_fint.o dp_tlong.o dp_flong.o sp_frexp.o sp_modf.o \
sp_div.o sp_mul.o sp_sub.o sp_add.o sp_fdp.o sp_cmp.o sp_logb.o \
sp_scalb.o sp_simple.o sp_tint.o sp_fint.o sp_tlong.o sp_flong.o \
dp_sqrt.o sp_sqrt.o kernel_linkage.o dsemul.o

1322
arch/mips/math-emu/cp1emu.c Normal file
View File

File diff suppressed because it is too large Load Diff

183
arch/mips/math-emu/dp_add.c Normal file
View File

@@ -0,0 +1,183 @@
/* IEEE754 floating point arithmetic
* double precision: common utilities
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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;
/* Infinity handling
*/
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 handling
*/
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);
}

67
arch/mips/math-emu/dp_cmp.c Normal file
View File

@@ -0,0 +1,67 @@
/* IEEE754 floating point arithmetic
* double precision: common utilities
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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, int sig)
{
COMPXDP;
COMPYDP;
EXPLODEXDP;
EXPLODEYDP;
FLUSHXDP;
FLUSHYDP;
CLEARCX; /* Even clear inexact flag here */
if (ieee754dp_isnan(x) || ieee754dp_isnan(y)) {
if (sig || 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 (sig && SETANDTESTCX(IEEE754_INVALID_OPERATION))
return ieee754si_xcpt(0, "fcmpf", x);
}
return 0;
} else {
s64 vx = x.bits;
s64 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;
}
}

157
arch/mips/math-emu/dp_div.c Normal file
View File

@@ -0,0 +1,157 @@
/* IEEE754 floating point arithmetic
* double precision: common utilities
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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 handling
*/
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 handling
*/
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 */
u64 rm = 0;
int re = xe - ye;
u64 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);
}
}

80
arch/mips/math-emu/dp_fint.c Normal file
View File

@@ -0,0 +1,80 @@
/* IEEE754 floating point arithmetic
* double precision: common utilities
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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_fint(int x)
{
COMPXDP;
CLEARCX;
xc = ( 0 ? xc : xc );
if (x == 0)
return ieee754dp_zero(0);
if (x == 1 || x == -1)
return ieee754dp_one(x < 0);
if (x == 10 || x == -10)
return ieee754dp_ten(x < 0);
xs = (x < 0);
if (xs) {
if (x == (1 << 31))
xm = ((unsigned) 1 << 31); /* max neg can't be safely negated */
else
xm = -x;
} else {
xm = x;
}
#if 1
/* normalize - result can never be inexact or overflow */
xe = DP_MBITS;
while ((xm >> DP_MBITS) == 0) {
xm <<= 1;
xe--;
}
return builddp(xs, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
#else
/* normalize */
xe = DP_MBITS + 3;
while ((xm >> (DP_MBITS + 3)) == 0) {
xm <<= 1;
xe--;
}
DPNORMRET1(xs, xe, xm, "fint", x);
#endif
}
ieee754dp ieee754dp_funs(unsigned int u)
{
if ((int) u < 0)
return ieee754dp_add(ieee754dp_1e31(),
ieee754dp_fint(u & ~(1 << 31)));
return ieee754dp_fint(u);
}

78
arch/mips/math-emu/dp_flong.c Normal file
View File

@@ -0,0 +1,78 @@
/* IEEE754 floating point arithmetic
* double precision: common utilities
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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_flong(s64 x)
{
COMPXDP;
CLEARCX;
xc = ( 0 ? xc : xc );
if (x == 0)
return ieee754dp_zero(0);
if (x == 1 || x == -1)
return ieee754dp_one(x < 0);
if (x == 10 || x == -10)
return ieee754dp_ten(x < 0);
xs = (x < 0);
if (xs) {
if (x == (1ULL << 63))
xm = (1ULL << 63); /* max neg can't be safely negated */
else
xm = -x;
} else {
xm = x;
}
/* normalize */
xe = DP_MBITS + 3;
if (xm >> (DP_MBITS + 1 + 3)) {
/* shunt out overflow bits */
while (xm >> (DP_MBITS + 1 + 3)) {
XDPSRSX1();
}
} else {
/* normalize in grs extended double precision */
while ((xm >> (DP_MBITS + 3)) == 0) {
xm <<= 1;
xe--;
}
}
DPNORMRET1(xs, xe, xm, "dp_flong", x);
}
ieee754dp ieee754dp_fulong(u64 u)
{
if ((s64) u < 0)
return ieee754dp_add(ieee754dp_1e63(),
ieee754dp_flong(u & ~(1ULL << 63)));
return ieee754dp_flong(u);
}

53
arch/mips/math-emu/dp_frexp.c Normal file
View File

@@ -0,0 +1,53 @@
/* IEEE754 floating point arithmetic
* double precision: common utilities
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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"
/* close to ieeep754dp_logb
*/
ieee754dp ieee754dp_frexp(ieee754dp x, int *eptr)
{
COMPXDP;
CLEARCX;
EXPLODEXDP;
switch (xc) {
case IEEE754_CLASS_SNAN:
case IEEE754_CLASS_QNAN:
case IEEE754_CLASS_INF:
case IEEE754_CLASS_ZERO:
*eptr = 0;
return x;
case IEEE754_CLASS_DNORM:
DPDNORMX;
break;
case IEEE754_CLASS_NORM:
break;
}
*eptr = xe + 1;
return builddp(xs, -1 + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
}

74
arch/mips/math-emu/dp_fsp.c Normal file
View File

@@ -0,0 +1,74 @@
/* IEEE754 floating point arithmetic
* double precision: common utilities
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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,
((u64) 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,
(u64) xm << (DP_MBITS - SP_MBITS));
}

54
arch/mips/math-emu/dp_logb.c Normal file
View File

@@ -0,0 +1,54 @@
/* IEEE754 floating point arithmetic
* double precision: common utilities
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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_logb(ieee754dp x)
{
COMPXDP;
CLEARCX;
EXPLODEXDP;
switch (xc) {
case IEEE754_CLASS_SNAN:
return ieee754dp_nanxcpt(x, "logb", x);
case IEEE754_CLASS_QNAN:
return x;
case IEEE754_CLASS_INF:
return ieee754dp_inf(0);
case IEEE754_CLASS_ZERO:
return ieee754dp_inf(1);
case IEEE754_CLASS_DNORM:
DPDNORMX;
break;
case IEEE754_CLASS_NORM:
break;
}
return ieee754dp_fint(xe);
}

80
arch/mips/math-emu/dp_modf.c Normal file
View File

@@ -0,0 +1,80 @@
/* IEEE754 floating point arithmetic
* double precision: common utilities
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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"
/* modf function is always exact for a finite number
*/
ieee754dp ieee754dp_modf(ieee754dp x, ieee754dp * ip)
{
COMPXDP;
CLEARCX;
EXPLODEXDP;
switch (xc) {
case IEEE754_CLASS_SNAN:
case IEEE754_CLASS_QNAN:
case IEEE754_CLASS_INF:
case IEEE754_CLASS_ZERO:
*ip = x;
return x;
case IEEE754_CLASS_DNORM:
/* far to small */
*ip = ieee754dp_zero(xs);
return x;
case IEEE754_CLASS_NORM:
break;
}
if (xe < 0) {
*ip = ieee754dp_zero(xs);
return x;
}
if (xe >= DP_MBITS) {
*ip = x;
return ieee754dp_zero(xs);
}
/* generate ipart mantissa by clearing bottom bits
*/
*ip = builddp(xs, xe + DP_EBIAS,
((xm >> (DP_MBITS - xe)) << (DP_MBITS - xe)) &
~DP_HIDDEN_BIT);
/* generate fpart mantissa by clearing top bits
* and normalizing (must be able to normalize)
*/
xm = (xm << (64 - (DP_MBITS - xe))) >> (64 - (DP_MBITS - xe));
if (xm == 0)
return ieee754dp_zero(xs);
while ((xm >> DP_MBITS) == 0) {
xm <<= 1;
xe--;
}
return builddp(xs, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
}

177
arch/mips/math-emu/dp_mul.c Normal file
View File

@@ -0,0 +1,177 @@
/* IEEE754 floating point arithmetic
* double precision: common utilities
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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 handling */
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;
u64 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) ((u64)(x) * (u64)y)
{
unsigned lxm = xm;
unsigned hxm = xm >> 32;
unsigned lym = ym;
unsigned hym = ym >> 32;
u64 lrm;
u64 hrm;
lrm = DPXMULT(lxm, lym);
hrm = DPXMULT(hxm, hym);
{
u64 t = DPXMULT(lxm, hym);
{
u64 at =
lrm + (t << 32);
hrm += at < lrm;
lrm = at;
}
hrm = hrm + (t >> 32);
}
{
u64 t = DPXMULT(hxm, lym);
{
u64 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 ((s64) 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);
}
}

58
arch/mips/math-emu/dp_scalb.c Normal file
View File

@@ -0,0 +1,58 @@
/* IEEE754 floating point arithmetic
* double precision: common utilities
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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_scalb(ieee754dp x, int n)
{
COMPXDP;
CLEARCX;
EXPLODEXDP;
switch (xc) {
case IEEE754_CLASS_SNAN:
return ieee754dp_nanxcpt(x, "scalb", x, n);
case IEEE754_CLASS_QNAN:
case IEEE754_CLASS_INF:
case IEEE754_CLASS_ZERO:
return x;
case IEEE754_CLASS_DNORM:
DPDNORMX;
break;
case IEEE754_CLASS_NORM:
break;
}
DPNORMRET2(xs, xe + n, xm << 3, "scalb", x, n);
}
ieee754dp ieee754dp_ldexp(ieee754dp x, int n)
{
return ieee754dp_scalb(x, n);
}

84
arch/mips/math-emu/dp_simple.c Normal file
View File

@@ -0,0 +1,84 @@
/* IEEE754 floating point arithmetic
* double precision: common utilities
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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;
}

165
arch/mips/math-emu/dp_sqrt.c Normal file
View File

@@ -0,0 +1,165 @@
/* IEEE754 floating point arithmetic
* double precision square root
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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 = ((u64) 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;
}

191
arch/mips/math-emu/dp_sub.c Normal file
View File

@@ -0,0 +1,191 @@
/* IEEE754 floating point arithmetic
* double precision: common utilities
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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_sub(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(), "sub", 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 handling
*/
case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
if (xs != ys)
return x;
SETCX(IEEE754_INVALID_OPERATION);
return ieee754dp_xcpt(ieee754dp_indef(), "sub", x, y);
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 ==
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):
/* quick fix up */
DPSIGN(y) ^= 1;
return y;
case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
DPDNORMX;
/* FAAL THOROUGH */
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
*/
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
*/
xm = xm + ym;
xe = xe;
xs = xs;
if (xm >> (DP_MBITS + 1 + 3)) { /* carry out */
xm = XDPSRS1(xm); /* shift preserving sticky */
xe++;
}
} else {
if (xm >= ym) {
xm = xm - ym;
xe = xe;
xs = xs;
} else {
xm = ym - xm;
xe = xe;
xs = ys;
}
if (xm == 0) {
if (ieee754_csr.rm == IEEE754_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_MBITS + 3)) == 0) {
xm <<= 1;
xe--;
}
}
DPNORMRET2(xs, xe, xm, "sub", x, y);
}

124
arch/mips/math-emu/dp_tint.c Normal file
View File

@@ -0,0 +1,124 @@
/* IEEE754 floating point arithmetic
* double precision: common utilities
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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 <linux/kernel.h>
#include "ieee754dp.h"
int ieee754dp_tint(ieee754dp x)
{
COMPXDP;
CLEARCX;
EXPLODEXDP;
FLUSHXDP;
switch (xc) {
case IEEE754_CLASS_SNAN:
case IEEE754_CLASS_QNAN:
case IEEE754_CLASS_INF:
SETCX(IEEE754_INVALID_OPERATION);
return ieee754si_xcpt(ieee754si_indef(), "dp_tint", x);
case IEEE754_CLASS_ZERO:
return 0;
case IEEE754_CLASS_DNORM:
case IEEE754_CLASS_NORM:
break;
}
if (xe > 31) {
/* Set invalid. We will only use overflow for floating
point overflow */
SETCX(IEEE754_INVALID_OPERATION);
return ieee754si_xcpt(ieee754si_indef(), "dp_tint", x);
}
/* oh gawd */
if (xe > DP_MBITS) {
xm <<= xe - DP_MBITS;
} else if (xe < DP_MBITS) {
u64 residue;
int round;
int sticky;
int odd;
if (xe < -1) {
residue = xm;
round = 0;
sticky = residue != 0;
xm = 0;
}
else {
residue = xm << (64 - DP_MBITS + xe);
round = (residue >> 63) != 0;
sticky = (residue << 1) != 0;
xm >>= DP_MBITS - xe;
}
/* Note: At this point upper 32 bits of xm are guaranteed
to be zero */
odd = (xm & 0x1) != 0x0;
switch (ieee754_csr.rm) {
case IEEE754_RN:
if (round && (sticky || odd))
xm++;
break;
case IEEE754_RZ:
break;
case IEEE754_RU: /* toward +Infinity */
if ((round || sticky) && !xs)
xm++;
break;
case IEEE754_RD: /* toward -Infinity */
if ((round || sticky) && xs)
xm++;
break;
}
/* look for valid corner case 0x80000000 */
if ((xm >> 31) != 0 && (xs == 0 || xm != 0x80000000)) {
/* This can happen after rounding */
SETCX(IEEE754_INVALID_OPERATION);
return ieee754si_xcpt(ieee754si_indef(), "dp_tint", x);
}
if (round || sticky)
SETCX(IEEE754_INEXACT);
}
if (xs)
return -xm;
else
return xm;
}
unsigned int ieee754dp_tuns(ieee754dp x)
{
ieee754dp hb = ieee754dp_1e31();
/* what if x < 0 ?? */
if (ieee754dp_lt(x, hb))
return (unsigned) ieee754dp_tint(x);
return (unsigned) ieee754dp_tint(ieee754dp_sub(x, hb)) |
((unsigned) 1 << 31);
}

127
arch/mips/math-emu/dp_tlong.c Normal file
View File

@@ -0,0 +1,127 @@
/* IEEE754 floating point arithmetic
* double precision: common utilities
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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"
s64 ieee754dp_tlong(ieee754dp x)
{
COMPXDP;
CLEARCX;
EXPLODEXDP;
FLUSHXDP;
switch (xc) {
case IEEE754_CLASS_SNAN:
case IEEE754_CLASS_QNAN:
case IEEE754_CLASS_INF:
SETCX(IEEE754_INVALID_OPERATION);
return ieee754di_xcpt(ieee754di_indef(), "dp_tlong", x);
case IEEE754_CLASS_ZERO:
return 0;
case IEEE754_CLASS_DNORM:
case IEEE754_CLASS_NORM:
break;
}
if (xe >= 63) {
/* look for valid corner case */
if (xe == 63 && xs && xm == DP_HIDDEN_BIT)
return -0x8000000000000000LL;
/* Set invalid. We will only use overflow for floating
point overflow */
SETCX(IEEE754_INVALID_OPERATION);
return ieee754di_xcpt(ieee754di_indef(), "dp_tlong", x);
}
/* oh gawd */
if (xe > DP_MBITS) {
xm <<= xe - DP_MBITS;
} else if (xe < DP_MBITS) {
u64 residue;
int round;
int sticky;
int odd;
if (xe < -1) {
residue = xm;
round = 0;
sticky = residue != 0;
xm = 0;
}
else {
/* Shifting a u64 64 times does not work,
* so we do it in two steps. Be aware that xe
* may be -1 */
residue = xm << (xe + 1);
residue <<= 63 - DP_MBITS;
round = (residue >> 63) != 0;
sticky = (residue << 1) != 0;
xm >>= DP_MBITS - xe;
}
odd = (xm & 0x1) != 0x0;
switch (ieee754_csr.rm) {
case IEEE754_RN:
if (round && (sticky || odd))
xm++;
break;
case IEEE754_RZ:
break;
case IEEE754_RU: /* toward +Infinity */
if ((round || sticky) && !xs)
xm++;
break;
case IEEE754_RD: /* toward -Infinity */
if ((round || sticky) && xs)
xm++;
break;
}
if ((xm >> 63) != 0) {
/* This can happen after rounding */
SETCX(IEEE754_INVALID_OPERATION);
return ieee754di_xcpt(ieee754di_indef(), "dp_tlong", x);
}
if (round || sticky)
SETCX(IEEE754_INEXACT);
}
if (xs)
return -xm;
else
return xm;
}
u64 ieee754dp_tulong(ieee754dp x)
{
ieee754dp hb = ieee754dp_1e63();
/* what if x < 0 ?? */
if (ieee754dp_lt(x, hb))
return (u64) ieee754dp_tlong(x);
return (u64) ieee754dp_tlong(ieee754dp_sub(x, hb)) |
(1ULL << 63);
}

172
arch/mips/math-emu/dsemul.c Normal file
View File

@@ -0,0 +1,172 @@
#include <linux/compiler.h>
#include <linux/mm.h>
#include <linux/signal.h>
#include <linux/smp.h>
#include <linux/smp_lock.h>
#include <asm/asm.h>
#include <asm/bootinfo.h>
#include <asm/byteorder.h>
#include <asm/cpu.h>
#include <asm/inst.h>
#include <asm/processor.h>
#include <asm/uaccess.h>
#include <asm/branch.h>
#include <asm/mipsregs.h>
#include <asm/system.h>
#include <asm/cacheflush.h>
#include <asm/fpu_emulator.h>
#include "ieee754.h"
#include "dsemul.h"
/* Strap kernel emulator for full MIPS IV emulation */
#ifdef __mips
#undef __mips
#endif
#define __mips 4
extern struct mips_fpu_emulator_private fpuemuprivate;
/*
* Emulate the arbritrary instruction ir at xcp->cp0_epc. Required when
* we have to emulate the instruction in a COP1 branch delay slot. Do
* not change cp0_epc due to the instruction
*
* According to the spec:
* 1) it shouldnt be a branch :-)
* 2) it can be a COP instruction :-(
* 3) if we are tring to run a protected memory space we must take
* special care on memory access instructions :-(
*/
/*
* "Trampoline" return routine to catch exception following
* execution of delay-slot instruction execution.
*/
struct emuframe {
mips_instruction emul;
mips_instruction badinst;
mips_instruction cookie;
gpreg_t epc;
};
int mips_dsemul(struct pt_regs *regs, mips_instruction ir, gpreg_t cpc)
{
extern asmlinkage void handle_dsemulret(void);
mips_instruction *dsemul_insns;
struct emuframe *fr;
int err;
if (ir == 0) { /* a nop is easy */
regs->cp0_epc = cpc;
regs->cp0_cause &= ~CAUSEF_BD;
return 0;
}
#ifdef DSEMUL_TRACE
printk("dsemul %lx %lx\n", regs->cp0_epc, cpc);
#endif
/*
* The strategy is to push the instruction onto the user stack
* and put a trap after it which we can catch and jump to
* the required address any alternative apart from full
* instruction emulation!!.
*
* Algorithmics used a system call instruction, and
* borrowed that vector. MIPS/Linux version is a bit
* more heavyweight in the interests of portability and
* multiprocessor support. For Linux we generate a
* an unaligned access and force an address error exception.
*
* For embedded systems (stand-alone) we prefer to use a
* non-existing CP1 instruction. This prevents us from emulating
* branches, but gives us a cleaner interface to the exception
* handler (single entry point).
*/
/* Ensure that the two instructions are in the same cache line */
dsemul_insns = (mips_instruction *) REG_TO_VA ((regs->regs[29] - sizeof(struct emuframe)) & ~0x7);
fr = (struct emuframe *) dsemul_insns;
/* Verify that the stack pointer is not competely insane */
if (unlikely(!access_ok(VERIFY_WRITE, fr, sizeof(struct emuframe))))
return SIGBUS;
err = __put_user(ir, &fr->emul);
err |= __put_user((mips_instruction)BADINST, &fr->badinst);
err |= __put_user((mips_instruction)BD_COOKIE, &fr->cookie);
err |= __put_user(cpc, &fr->epc);
if (unlikely(err)) {
fpuemuprivate.stats.errors++;
return SIGBUS;
}
regs->cp0_epc = VA_TO_REG & fr->emul;
flush_cache_sigtramp((unsigned long)&fr->badinst);
return SIGILL; /* force out of emulation loop */
}
int do_dsemulret(struct pt_regs *xcp)
{
struct emuframe *fr;
gpreg_t epc;
u32 insn, cookie;
int err = 0;
fr = (struct emuframe *) (xcp->cp0_epc - sizeof(mips_instruction));
/*
* If we can't even access the area, something is very wrong, but we'll
* leave that to the default handling
*/
if (!access_ok(VERIFY_READ, fr, sizeof(struct emuframe)))
return 0;
/*
* Do some sanity checking on the stackframe:
*
* - Is the instruction pointed to by the EPC an BADINST?
* - Is the following memory word the BD_COOKIE?
*/
err = __get_user(insn, &fr->badinst);
err |= __get_user(cookie, &fr->cookie);
if (unlikely(err || (insn != BADINST) || (cookie != BD_COOKIE))) {
fpuemuprivate.stats.errors++;
return 0;
}
/*
* At this point, we are satisfied that it's a BD emulation trap. Yes,
* a user might have deliberately put two malformed and useless
* instructions in a row in his program, in which case he's in for a
* nasty surprise - the next instruction will be treated as a
* continuation address! Alas, this seems to be the only way that we
* can handle signals, recursion, and longjmps() in the context of
* emulating the branch delay instruction.
*/
#ifdef DSEMUL_TRACE
printk("dsemulret\n");
#endif
if (__get_user(epc, &fr->epc)) { /* Saved EPC */
/* This is not a good situation to be in */
force_sig(SIGBUS, current);
return 0;
}
/* Set EPC to return to post-branch instruction */
xcp->cp0_epc = epc;
return 1;
}

23
arch/mips/math-emu/dsemul.h Normal file
View File

@@ -0,0 +1,23 @@
typedef long gpreg_t;
typedef void *vaddr_t;
#define REG_TO_VA (vaddr_t)
#define VA_TO_REG (gpreg_t)
int mips_dsemul(struct pt_regs *regs, mips_instruction ir, gpreg_t cpc);
int do_dsemulret(struct pt_regs *xcp);
/* Instruction which will always cause an address error */
#define AdELOAD 0x8c000001 /* lw $0,1($0) */
/* Instruction which will plainly cause a CP1 exception when FPU is disabled */
#define CP1UNDEF 0x44400001 /* cfc1 $0,$0 undef */
/* Instruction inserted following the badinst to further tag the sequence */
#define BD_COOKIE 0x0000bd36 /* tne $0,$0 with baggage */
/* Setup which instruction to use for trampoline */
#ifdef STANDALONE_EMULATOR
#define BADINST CP1UNDEF
#else
#define BADINST AdELOAD
#endif

138
arch/mips/math-emu/ieee754.c Normal file
View File

@@ -0,0 +1,138 @@
/* ieee754 floating point arithmetic
* single and double precision
*
* BUGS
* not much dp done
* doesn't generate IEEE754_INEXACT
*
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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;
}
s64 ieee754di_xcpt(s64 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;
}

489
arch/mips/math-emu/ieee754.h Normal file
View File

@@ -0,0 +1,489 @@
/* single and double precision fp ops
* missing extended precision.
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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, kevink@mips.com and Carsten Langgard, carstenl@mips.com
* 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 {
u64 mant:52;
unsigned int bexp:11;
unsigned int sign:1;
} parts;
u64 bits;
double d;
} ieee754dp;
typedef union _ieee754sp {
struct ieee754sp_konst parts;
float f;
u32 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;
u64 mant:52;
} parts;
double d;
u64 bits;
} ieee754dp;
struct ieee754sp_konst {
unsigned sign:1;
unsigned bexp:8;
unsigned mant:23;
};
typedef union _ieee754sp {
struct ieee754sp_konst parts;
float f;
u32 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(s64 x);
ieee754sp ieee754sp_fulong(u64 x);
ieee754sp ieee754sp_fdp(ieee754dp x);
int ieee754sp_tint(ieee754sp x);
unsigned int ieee754sp_tuns(ieee754sp x);
s64 ieee754sp_tlong(ieee754sp x);
u64 ieee754sp_tulong(ieee754sp x);
int ieee754sp_cmp(ieee754sp x, ieee754sp y, int cop, int sig);
/*
* 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(s64 x);
ieee754dp ieee754dp_fulong(u64 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);
s64 ieee754dp_tlong(ieee754dp x);
u64 ieee754dp_tulong(ieee754dp x);
int ieee754dp_cmp(ieee754dp x, ieee754dp y, int cop, int sig);
/*
* 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, 0);
}
static __inline int ieee754sp_ne(ieee754sp x, ieee754sp y)
{
return ieee754sp_cmp(x, y,
IEEE754_CLT | IEEE754_CGT | IEEE754_CUN, 0);
}
static __inline int ieee754sp_lt(ieee754sp x, ieee754sp y)
{
return ieee754sp_cmp(x, y, IEEE754_CLT, 0);
}
static __inline int ieee754sp_le(ieee754sp x, ieee754sp y)
{
return ieee754sp_cmp(x, y, IEEE754_CLT | IEEE754_CEQ, 0);
}
static __inline int ieee754sp_gt(ieee754sp x, ieee754sp y)
{
return ieee754sp_cmp(x, y, IEEE754_CGT, 0);
}
static __inline int ieee754sp_ge(ieee754sp x, ieee754sp y)
{
return ieee754sp_cmp(x, y, IEEE754_CGT | IEEE754_CEQ, 0);
}
static __inline int ieee754dp_eq(ieee754dp x, ieee754dp y)
{
return ieee754dp_cmp(x, y, IEEE754_CEQ, 0);
}
static __inline int ieee754dp_ne(ieee754dp x, ieee754dp y)
{
return ieee754dp_cmp(x, y,
IEEE754_CLT | IEEE754_CGT | IEEE754_CUN, 0);
}
static __inline int ieee754dp_lt(ieee754dp x, ieee754dp y)
{
return ieee754dp_cmp(x, y, IEEE754_CLT, 0);
}
static __inline int ieee754dp_le(ieee754dp x, ieee754dp y)
{
return ieee754dp_cmp(x, y, IEEE754_CLT | IEEE754_CEQ, 0);
}
static __inline int ieee754dp_gt(ieee754dp x, ieee754dp y)
{
return ieee754dp_cmp(x, y, IEEE754_CGT, 0);
}
static __inline int ieee754dp_ge(ieee754dp x, ieee754dp y)
{
return ieee754dp_cmp(x, y, IEEE754_CGT | IEEE754_CEQ, 0);
}
/* 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() ((s64)(~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) */
s64 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)

138
arch/mips/math-emu/ieee754d.c Normal file
View File

@@ -0,0 +1,138 @@
/*
* Some debug functions
*
* MIPS floating point support
*
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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
* Modified to build and operate in Linux kernel environment.
*
* Kevin D. Kissell, kevink@mips.com and Carsten Langgaard, carstenl@mips.com
* Copyright (C) 2000 MIPS Technologies, Inc. All rights reserved.
*/
#include <linux/kernel.h>
#include "ieee754.h"
#define DP_EBIAS 1023
#define DP_EMIN (-1022)
#define DP_EMAX 1023
#define DP_FBITS 52
#define SP_EBIAS 127
#define SP_EMIN (-126)
#define SP_EMAX 127
#define SP_FBITS 23
#define DP_MBIT(x) ((u64)1 << (x))
#define DP_HIDDEN_BIT DP_MBIT(DP_FBITS)
#define DP_SIGN_BIT DP_MBIT(63)
#define SP_MBIT(x) ((u32)1 << (x))
#define SP_HIDDEN_BIT SP_MBIT(SP_FBITS)
#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)
ieee754dp ieee754dp_dump(char *m, ieee754dp x)
{
int i;
printk("%s", m);
printk("<%08x,%08x>\n", (unsigned) (x.bits >> 32),
(unsigned) x.bits);
printk("\t=");
switch (ieee754dp_class(x)) {
case IEEE754_CLASS_QNAN:
case IEEE754_CLASS_SNAN:
printk("Nan %c", DPSIGN(x) ? '-' : '+');
for (i = DP_FBITS - 1; i >= 0; i--)
printk("%c", DPMANT(x) & DP_MBIT(i) ? '1' : '0');
break;
case IEEE754_CLASS_INF:
printk("%cInfinity", DPSIGN(x) ? '-' : '+');
break;
case IEEE754_CLASS_ZERO:
printk("%cZero", DPSIGN(x) ? '-' : '+');
break;
case IEEE754_CLASS_DNORM:
printk("%c0.", DPSIGN(x) ? '-' : '+');
for (i = DP_FBITS - 1; i >= 0; i--)
printk("%c", DPMANT(x) & DP_MBIT(i) ? '1' : '0');
printk("e%d", DPBEXP(x) - DP_EBIAS);
break;
case IEEE754_CLASS_NORM:
printk("%c1.", DPSIGN(x) ? '-' : '+');
for (i = DP_FBITS - 1; i >= 0; i--)
printk("%c", DPMANT(x) & DP_MBIT(i) ? '1' : '0');
printk("e%d", DPBEXP(x) - DP_EBIAS);
break;
default:
printk("Illegal/Unknown IEEE754 value class");
}
printk("\n");
return x;
}
ieee754sp ieee754sp_dump(char *m, ieee754sp x)
{
int i;
printk("%s=", m);
printk("<%08x>\n", (unsigned) x.bits);
printk("\t=");
switch (ieee754sp_class(x)) {
case IEEE754_CLASS_QNAN:
case IEEE754_CLASS_SNAN:
printk("Nan %c", SPSIGN(x) ? '-' : '+');
for (i = SP_FBITS - 1; i >= 0; i--)
printk("%c", SPMANT(x) & SP_MBIT(i) ? '1' : '0');
break;
case IEEE754_CLASS_INF:
printk("%cInfinity", SPSIGN(x) ? '-' : '+');
break;
case IEEE754_CLASS_ZERO:
printk("%cZero", SPSIGN(x) ? '-' : '+');
break;
case IEEE754_CLASS_DNORM:
printk("%c0.", SPSIGN(x) ? '-' : '+');
for (i = SP_FBITS - 1; i >= 0; i--)
printk("%c", SPMANT(x) & SP_MBIT(i) ? '1' : '0');
printk("e%d", SPBEXP(x) - SP_EBIAS);
break;
case IEEE754_CLASS_NORM:
printk("%c1.", SPSIGN(x) ? '-' : '+');
for (i = SP_FBITS - 1; i >= 0; i--)
printk("%c", SPMANT(x) & SP_MBIT(i) ? '1' : '0');
printk("e%d", SPBEXP(x) - SP_EBIAS);
break;
default:
printk("Illegal/Unknown IEEE754 value class");
}
printk("\n");
return x;
}

243
arch/mips/math-emu/ieee754dp.c Normal file
View File

@@ -0,0 +1,243 @@
/* IEEE754 floating point arithmetic
* double precision: common utilities
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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 (!SETANDTESTCX(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;
}
static u64 get_rounding(int sn, u64 xm)
{
/* inexact must round of 3 bits
*/
if (xm & (DP_MBIT(3) - 1)) {
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;
}
}
return xm;
}
/* generate a normal/denormal number with over,under handling
* sn is sign
* xe is an unbiased exponent
* xm is 3bit extended precision value.
*/
ieee754dp ieee754dp_format(int sn, int xe, u64 xm)
{
assert(xm); /* we don't 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);
}
}
if (xe == DP_EMIN - 1
&& get_rounding(sn, xm) >> (DP_MBITS + 1 + 3))
{
/* Not tiny after rounding */
SETCX(IEEE754_INEXACT);
xm = get_rounding(sn, xm);
xm >>= 1;
/* Clear grs bits */
xm &= ~(DP_MBIT(3) - 1);
xe++;
}
else {
/* 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);
if ((xm & (DP_HIDDEN_BIT << 3)) == 0) {
SETCX(IEEE754_UNDERFLOW);
}
/* inexact must round of 3 bits
*/
xm = get_rounding(sn, xm);
/* 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);
if (ieee754_csr.mx & IEEE754_UNDERFLOW)
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);
}
}

83
arch/mips/math-emu/ieee754dp.h Normal file
View File

@@ -0,0 +1,83 @@
/*
* IEEE754 floating point
* double precision internal header file
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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 assert(expr) ((void)0)
/* 3bit extended double precision sticky right shift */
#define XDPSRS(v,rs) \
((rs > (DP_MBITS+3))?1:((v) >> (rs)) | ((v) << (64-(rs)) != 0))
#define XDPSRSX1() \
(xe++, (xm = (xm >> 1) | (xm & 1)))
#define XDPSRS1(v) \
(((v) >> 1) | ((v) & 1))
/* convert denormal to normalized with extended exponent */
#define DPDNORMx(m,e) \
while( (m >> DP_MBITS) == 0) { m <<= 1; e--; }
#define DPDNORMX DPDNORMx(xm,xe)
#define DPDNORMY DPDNORMx(ym,ye)
static __inline ieee754dp builddp(int s, int bx, u64 m)
{
ieee754dp r;
assert((s) == 0 || (s) == 1);
assert((bx) >= DP_EMIN - 1 + DP_EBIAS
&& (bx) <= DP_EMAX + 1 + DP_EBIAS);
assert(((m) >> DP_MBITS) == 0);
r.parts.sign = s;
r.parts.bexp = bx;
r.parts.mant = m;
return r;
}
extern int ieee754dp_isnan(ieee754dp);
extern int ieee754dp_issnan(ieee754dp);
extern int ieee754si_xcpt(int, const char *, ...);
extern s64 ieee754di_xcpt(s64, const char *, ...);
extern ieee754dp ieee754dp_xcpt(ieee754dp, const char *, ...);
extern ieee754dp ieee754dp_nanxcpt(ieee754dp, const char *, ...);
extern ieee754dp ieee754dp_bestnan(ieee754dp, ieee754dp);
extern ieee754dp ieee754dp_format(int, int, u64);
#define DPNORMRET2(s,e,m,name,a0,a1) \
{ \
ieee754dp V = ieee754dp_format(s,e,m); \
if(TSTX()) \
return ieee754dp_xcpt(V,name,a0,a1); \
else \
return V; \
}
#define DPNORMRET1(s,e,m,name,a0) DPNORMRET2(s,e,m,name,a0,a0)

165
arch/mips/math-emu/ieee754int.h Normal file
View File

@@ -0,0 +1,165 @@
/*
* IEEE754 floating point
* common internal header file
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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) ((u64)1 << (x))
#define DP_HIDDEN_BIT DP_MBIT(DP_MBITS)
#define DP_SIGN_BIT DP_MBIT(63)
#define SP_MBIT(x) ((u32)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))
#define SETANDTESTCX(x) \
(SETCX(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 \
u64 xm; int xe; int xs; int xc
#define COMPYDP \
u64 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)

56
arch/mips/math-emu/ieee754m.c Normal file
View File

@@ -0,0 +1,56 @@
/*
* floor, trunc, ceil
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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"
ieee754dp ieee754dp_floor(ieee754dp x)
{
ieee754dp i;
if (ieee754dp_lt(ieee754dp_modf(x, &i), ieee754dp_zero(0)))
return ieee754dp_sub(i, ieee754dp_one(0));
else
return i;
}
ieee754dp ieee754dp_ceil(ieee754dp x)
{
ieee754dp i;
if (ieee754dp_gt(ieee754dp_modf(x, &i), ieee754dp_zero(0)))
return ieee754dp_add(i, ieee754dp_one(0));
else
return i;
}
ieee754dp ieee754dp_trunc(ieee754dp x)
{
ieee754dp i;
(void) ieee754dp_modf(x, &i);
return i;
}

243
arch/mips/math-emu/ieee754sp.c Normal file
View File

@@ -0,0 +1,243 @@
/* IEEE754 floating point arithmetic
* single precision
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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 (!SETANDTESTCX(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;
}
static unsigned get_rounding(int sn, unsigned xm)
{
/* inexact must round of 3 bits
*/
if (xm & (SP_MBIT(3) - 1)) {
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;
}
}
return xm;
}
/* generate a normal/denormal number with over,under handling
* 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 don't 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);
}
}
if (xe == SP_EMIN - 1
&& get_rounding(sn, xm) >> (SP_MBITS + 1 + 3))
{
/* Not tiny after rounding */
SETCX(IEEE754_INEXACT);
xm = get_rounding(sn, xm);
xm >>= 1;
/* Clear grs bits */
xm &= ~(SP_MBIT(3) - 1);
xe++;
}
else {
/* 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);
if ((xm & (SP_HIDDEN_BIT << 3)) == 0) {
SETCX(IEEE754_UNDERFLOW);
}
/* inexact must round of 3 bits
*/
xm = get_rounding(sn, xm);
/* 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);
if (ieee754_csr.mx & IEEE754_UNDERFLOW)
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);
}
}

89
arch/mips/math-emu/ieee754sp.h Normal file
View File

@@ -0,0 +1,89 @@
/*
* IEEE754 floating point
* double precision internal header file
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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 assert(expr) ((void)0)
/* 3bit extended single precision sticky right shift */
#define SPXSRSXn(rs) \
(xe += rs, \
xm = (rs > (SP_MBITS+3))?1:((xm) >> (rs)) | ((xm) << (32-(rs)) != 0))
#define SPXSRSX1() \
(xe++, (xm = (xm >> 1) | (xm & 1)))
#define SPXSRSYn(rs) \
(ye+=rs, \
ym = (rs > (SP_MBITS+3))?1:((ym) >> (rs)) | ((ym) << (32-(rs)) != 0))
#define SPXSRSY1() \
(ye++, (ym = (ym >> 1) | (ym & 1)))
/* convert denormal to normalized with extended exponent */
#define SPDNORMx(m,e) \
while( (m >> SP_MBITS) == 0) { m <<= 1; e--; }
#define SPDNORMX SPDNORMx(xm,xe)
#define SPDNORMY SPDNORMx(ym,ye)
static __inline ieee754sp buildsp(int s, int bx, unsigned m)
{
ieee754sp r;
assert((s) == 0 || (s) == 1);
assert((bx) >= SP_EMIN - 1 + SP_EBIAS
&& (bx) <= SP_EMAX + 1 + SP_EBIAS);
assert(((m) >> SP_MBITS) == 0);
r.parts.sign = s;
r.parts.bexp = bx;
r.parts.mant = m;
return r;
}
extern int ieee754sp_isnan(ieee754sp);
extern int ieee754sp_issnan(ieee754sp);
extern int ieee754si_xcpt(int, const char *, ...);
extern s64 ieee754di_xcpt(s64, const char *, ...);
extern ieee754sp ieee754sp_xcpt(ieee754sp, const char *, ...);
extern ieee754sp ieee754sp_nanxcpt(ieee754sp, const char *, ...);
extern ieee754sp ieee754sp_bestnan(ieee754sp, ieee754sp);
extern ieee754sp ieee754sp_format(int, int, unsigned);
#define SPNORMRET2(s,e,m,name,a0,a1) \
{ \
ieee754sp V = ieee754sp_format(s,e,m); \
if(TSTX()) \
return ieee754sp_xcpt(V,name,a0,a1); \
else \
return V; \
}
#define SPNORMRET1(s,e,m,name,a0) SPNORMRET2(s,e,m,name,a0,a0)

49
arch/mips/math-emu/ieee754xcpt.c Normal file
View File

@@ -0,0 +1,49 @@
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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
* Added preprocessor hacks to map to Linux kernel diagnostics.
*
* Kevin D. Kissell, kevink@mips.com and Carsten Langgaard, carstenl@mips.com
* Copyright (C) 2000 MIPS Technologies, Inc. All rights reserved.
*************************************************************************/
#include <linux/kernel.h>
#include "ieee754.h"
/*
* Very naff exception handler (you can plug in your own and
* override this).
*/
static const char *const rtnames[] = {
"sp", "dp", "xp", "si", "di"
};
void ieee754_xcpt(struct ieee754xctx *xcp)
{
printk(KERN_DEBUG "floating point exception in \"%s\", type=%s\n",
xcp->op, rtnames[xcp->rt]);
}

125
arch/mips/math-emu/kernel_linkage.c Normal file
View File

@@ -0,0 +1,125 @@
/*
* Kevin D. Kissell, kevink@mips and Carsten Langgaard, carstenl@mips.com
* 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/config.h>
#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.fcr31 = 0;
for (i = 0; i < 32; i++) {
current->thread.fpu.soft.fpr[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.fpr[i],
&sc->sc_fpregs[i]);
}
err |= __put_user(current->thread.fpu.soft.fcr31, &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.fpr[i],
&sc->sc_fpregs[i]);
}
err |= __get_user(current->thread.fpu.soft.fcr31, &sc->sc_fpc_csr);
err |= __get_user(fpuemuprivate.eir, &sc->sc_fpc_eir);
return err;
}
#ifdef CONFIG_MIPS64
/*
* This is the o32 version
*/
int fpu_emulator_save_context32(struct sigcontext32 *sc)
{
int i;
int err = 0;
for (i = 0; i < 32; i+=2) {
err |=
__put_user(current->thread.fpu.soft.fpr[i],
&sc->sc_fpregs[i]);
}
err |= __put_user(current->thread.fpu.soft.fcr31, &sc->sc_fpc_csr);
err |= __put_user(fpuemuprivate.eir, &sc->sc_fpc_eir);
return err;
}
int fpu_emulator_restore_context32(struct sigcontext32 *sc)
{
int i;
int err = 0;
for (i = 0; i < 32; i+=2) {
err |=
__get_user(current->thread.fpu.soft.fpr[i],
&sc->sc_fpregs[i]);
}
err |= __get_user(current->thread.fpu.soft.fcr31, &sc->sc_fpc_csr);
err |= __get_user(fpuemuprivate.eir, &sc->sc_fpc_eir);
return err;
}
#endif

177
arch/mips/math-emu/sp_add.c Normal file
View File

@@ -0,0 +1,177 @@
/* IEEE754 floating point arithmetic
* single precision
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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):
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;
SETCX(IEEE754_INVALID_OPERATION);
return ieee754sp_xcpt(ieee754sp_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 handling
*/
case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
if (xs == ys)
return x;
else
return ieee754sp_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):
SPDNORMX;
case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
SPDNORMY;
break;
case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
SPDNORMX;
break;
case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
break;
}
assert(xm & SP_HIDDEN_BIT);
assert(ym & SP_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;
SPXSRSYn(s);
} else if (ye > xe) {
/* have to shift x fraction right to align
*/
int s = ye - xe;
SPXSRSXn(s);
}
assert(xe == ye);
assert(xe <= SP_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 >> (SP_MBITS + 1 + 3)) { /* carry out */
SPXSRSX1();
}
} else {
if (xm >= ym) {
xm = xm - ym;
xe = xe;
xs = xs;
} else {
xm = ym - xm;
xe = xe;
xs = ys;
}
if (xm == 0)
return ieee754sp_zero(ieee754_csr.rm ==
IEEE754_RD);
/* normalize in extended single precision */
while ((xm >> (SP_MBITS + 3)) == 0) {
xm <<= 1;
xe--;
}
}
SPNORMRET2(xs, xe, xm, "add", x, y);
}

67
arch/mips/math-emu/sp_cmp.c Normal file
View File

@@ -0,0 +1,67 @@
/* IEEE754 floating point arithmetic
* single precision
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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_cmp(ieee754sp x, ieee754sp y, int cmp, int sig)
{
COMPXSP;
COMPYSP;
EXPLODEXSP;
EXPLODEYSP;
FLUSHXSP;
FLUSHYSP;
CLEARCX; /* Even clear inexact flag here */
if (ieee754sp_isnan(x) || ieee754sp_isnan(y)) {
if (sig || 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 (sig && SETANDTESTCX(IEEE754_INVALID_OPERATION))
return ieee754si_xcpt(0, "fcmpf", x);
}
return 0;
} else {
int vx = x.bits;
int vy = y.bits;
if (vx < 0)
vx = -vx ^ SP_SIGN_BIT;
if (vy < 0)
vy = -vy ^ SP_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;
}
}

157
arch/mips/math-emu/sp_div.c Normal file
View File

@@ -0,0 +1,157 @@
/* IEEE754 floating point arithmetic
* single precision
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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_div(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(), "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 handling
*/
case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
SETCX(IEEE754_INVALID_OPERATION);
return ieee754sp_xcpt(ieee754sp_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 ieee754sp_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 ieee754sp_inf(xs ^ ys);
/* Zero handling
*/
case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
SETCX(IEEE754_INVALID_OPERATION);
return ieee754sp_xcpt(ieee754sp_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 ieee754sp_xcpt(ieee754sp_inf(xs ^ ys), "div", x, y);
case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
return ieee754sp_zero(xs == ys ? 0 : 1);
case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
SPDNORMX;
case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
SPDNORMY;
break;
case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
SPDNORMX;
break;
case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
break;
}
assert(xm & SP_HIDDEN_BIT);
assert(ym & SP_HIDDEN_BIT);
/* provide rounding space */
xm <<= 3;
ym <<= 3;
{
/* now the dirty work */
unsigned rm = 0;
int re = xe - ye;
unsigned bm;
for (bm = SP_MBIT(SP_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 >> (SP_MBITS + 3)) == 0) {
rm <<= 1;
re--;
}
SPNORMRET2(xs == ys ? 0 : 1, re, rm, "div", x, y);
}
}

77
arch/mips/math-emu/sp_fdp.c Normal file
View File

@@ -0,0 +1,77 @@
/* IEEE754 floating point arithmetic
* single precision
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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_fdp(ieee754dp x)
{
COMPXDP;
ieee754sp nan;
EXPLODEXDP;
CLEARCX;
FLUSHXDP;
switch (xc) {
case IEEE754_CLASS_SNAN:
SETCX(IEEE754_INVALID_OPERATION);
return ieee754sp_nanxcpt(ieee754sp_indef(), "fdp");
case IEEE754_CLASS_QNAN:
nan = buildsp(xs, SP_EMAX + 1 + SP_EBIAS, (u32)
(xm >> (DP_MBITS - SP_MBITS)));
if (!ieee754sp_isnan(nan))
nan = ieee754sp_indef();
return ieee754sp_nanxcpt(nan, "fdp", x);
case IEEE754_CLASS_INF:
return ieee754sp_inf(xs);
case IEEE754_CLASS_ZERO:
return ieee754sp_zero(xs);
case IEEE754_CLASS_DNORM:
/* can't possibly be sp representable */
SETCX(IEEE754_UNDERFLOW);
SETCX(IEEE754_INEXACT);
if ((ieee754_csr.rm == IEEE754_RU && !xs) ||
(ieee754_csr.rm == IEEE754_RD && xs))
return ieee754sp_xcpt(ieee754sp_mind(xs), "fdp", x);
return ieee754sp_xcpt(ieee754sp_zero(xs), "fdp", x);
case IEEE754_CLASS_NORM:
break;
}
{
u32 rm;
/* convert from DP_MBITS to SP_MBITS+3 with sticky right shift
*/
rm = (xm >> (DP_MBITS - (SP_MBITS + 3))) |
((xm << (64 - (DP_MBITS - (SP_MBITS + 3)))) != 0);
SPNORMRET1(xs, xe, rm, "fdp", x);
}
}

80
arch/mips/math-emu/sp_fint.c Normal file
View File

@@ -0,0 +1,80 @@
/* IEEE754 floating point arithmetic
* single precision
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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_fint(int x)
{
COMPXSP;
CLEARCX;
xc = ( 0 ? xc : xc );
if (x == 0)
return ieee754sp_zero(0);
if (x == 1 || x == -1)
return ieee754sp_one(x < 0);
if (x == 10 || x == -10)
return ieee754sp_ten(x < 0);
xs = (x < 0);
if (xs) {
if (x == (1 << 31))
xm = ((unsigned) 1 << 31); /* max neg can't be safely negated */
else
xm = -x;
} else {
xm = x;
}
xe = SP_MBITS + 3;
if (xm >> (SP_MBITS + 1 + 3)) {
/* shunt out overflow bits
*/
while (xm >> (SP_MBITS + 1 + 3)) {
SPXSRSX1();
}
} else {
/* normalize in grs extended single precision
*/
while ((xm >> (SP_MBITS + 3)) == 0) {
xm <<= 1;
xe--;
}
}
SPNORMRET1(xs, xe, xm, "fint", x);
}
ieee754sp ieee754sp_funs(unsigned int u)
{
if ((int) u < 0)
return ieee754sp_add(ieee754sp_1e31(),
ieee754sp_fint(u & ~(1 << 31)));
return ieee754sp_fint(u);
}

79
arch/mips/math-emu/sp_flong.c Normal file
View File

@@ -0,0 +1,79 @@
/* IEEE754 floating point arithmetic
* single precision
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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_flong(s64 x)
{
COMPXDP; /* <--- need 64-bit mantissa temp */
CLEARCX;
xc = ( 0 ? xc : xc );
if (x == 0)
return ieee754sp_zero(0);
if (x == 1 || x == -1)
return ieee754sp_one(x < 0);
if (x == 10 || x == -10)
return ieee754sp_ten(x < 0);
xs = (x < 0);
if (xs) {
if (x == (1ULL << 63))
xm = (1ULL << 63); /* max neg can't be safely negated */
else
xm = -x;
} else {
xm = x;
}
xe = SP_MBITS + 3;
if (xm >> (SP_MBITS + 1 + 3)) {
/* shunt out overflow bits
*/
while (xm >> (SP_MBITS + 1 + 3)) {
SPXSRSX1();
}
} else {
/* normalize in grs extended single precision */
while ((xm >> (SP_MBITS + 3)) == 0) {
xm <<= 1;
xe--;
}
}
SPNORMRET1(xs, xe, xm, "sp_flong", x);
}
ieee754sp ieee754sp_fulong(u64 u)
{
if ((s64) u < 0)
return ieee754sp_add(ieee754sp_1e63(),
ieee754sp_flong(u & ~(1ULL << 63)));
return ieee754sp_flong(u);
}

53
arch/mips/math-emu/sp_frexp.c Normal file
View File

@@ -0,0 +1,53 @@
/* IEEE754 floating point arithmetic
* single precision
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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"
/* close to ieeep754sp_logb
*/
ieee754sp ieee754sp_frexp(ieee754sp x, int *eptr)
{
COMPXSP;
CLEARCX;
EXPLODEXSP;
switch (xc) {
case IEEE754_CLASS_SNAN:
case IEEE754_CLASS_QNAN:
case IEEE754_CLASS_INF:
case IEEE754_CLASS_ZERO:
*eptr = 0;
return x;
case IEEE754_CLASS_DNORM:
SPDNORMX;
break;
case IEEE754_CLASS_NORM:
break;
}
*eptr = xe + 1;
return buildsp(xs, -1 + SP_EBIAS, xm & ~SP_HIDDEN_BIT);
}

54
arch/mips/math-emu/sp_logb.c Normal file
View File

@@ -0,0 +1,54 @@
/* IEEE754 floating point arithmetic
* single precision
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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_logb(ieee754sp x)
{
COMPXSP;
CLEARCX;
EXPLODEXSP;
switch (xc) {
case IEEE754_CLASS_SNAN:
return ieee754sp_nanxcpt(x, "logb", x);
case IEEE754_CLASS_QNAN:
return x;
case IEEE754_CLASS_INF:
return ieee754sp_inf(0);
case IEEE754_CLASS_ZERO:
return ieee754sp_inf(1);
case IEEE754_CLASS_DNORM:
SPDNORMX;
break;
case IEEE754_CLASS_NORM:
break;
}
return ieee754sp_fint(xe);
}

80
arch/mips/math-emu/sp_modf.c Normal file
View File

@@ -0,0 +1,80 @@
/* IEEE754 floating point arithmetic
* single precision
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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"
/* modf function is always exact for a finite number
*/
ieee754sp ieee754sp_modf(ieee754sp x, ieee754sp * ip)
{
COMPXSP;
CLEARCX;
EXPLODEXSP;
switch (xc) {
case IEEE754_CLASS_SNAN:
case IEEE754_CLASS_QNAN:
case IEEE754_CLASS_INF:
case IEEE754_CLASS_ZERO:
*ip = x;
return x;
case IEEE754_CLASS_DNORM:
/* far to small */
*ip = ieee754sp_zero(xs);
return x;
case IEEE754_CLASS_NORM:
break;
}
if (xe < 0) {
*ip = ieee754sp_zero(xs);
return x;
}
if (xe >= SP_MBITS) {
*ip = x;
return ieee754sp_zero(xs);
}
/* generate ipart mantissa by clearing bottom bits
*/
*ip = buildsp(xs, xe + SP_EBIAS,
((xm >> (SP_MBITS - xe)) << (SP_MBITS - xe)) &
~SP_HIDDEN_BIT);
/* generate fpart mantissa by clearing top bits
* and normalizing (must be able to normalize)
*/
xm = (xm << (32 - (SP_MBITS - xe))) >> (32 - (SP_MBITS - xe));
if (xm == 0)
return ieee754sp_zero(xs);
while ((xm >> SP_MBITS) == 0) {
xm <<= 1;
xe--;
}
return buildsp(xs, xe + SP_EBIAS, xm & ~SP_HIDDEN_BIT);
}

171
arch/mips/math-emu/sp_mul.c Normal file
View File

@@ -0,0 +1,171 @@
/* IEEE754 floating point arithmetic
* single precision
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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_mul(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(), "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 handling */
case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
SETCX(IEEE754_INVALID_OPERATION);
return ieee754sp_xcpt(ieee754sp_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 ieee754sp_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 ieee754sp_zero(xs ^ ys);
case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
SPDNORMX;
case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
SPDNORMY;
break;
case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
SPDNORMX;
break;
case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
break;
}
/* rm = xm * ym, re = xe+ye basicly */
assert(xm & SP_HIDDEN_BIT);
assert(ym & SP_HIDDEN_BIT);
{
int re = xe + ye;
int rs = xs ^ ys;
unsigned rm;
/* shunt to top of word */
xm <<= 32 - (SP_MBITS + 1);
ym <<= 32 - (SP_MBITS + 1);
/* multiply 32bits xm,ym to give high 32bits rm with stickness
*/
{
unsigned short lxm = xm & 0xffff;
unsigned short hxm = xm >> 16;
unsigned short lym = ym & 0xffff;
unsigned short hym = ym >> 16;
unsigned lrm;
unsigned hrm;
lrm = lxm * lym; /* 16 * 16 => 32 */
hrm = hxm * hym; /* 16 * 16 => 32 */
{
unsigned t = lxm * hym; /* 16 * 16 => 32 */
{
unsigned at = lrm + (t << 16);
hrm += at < lrm;
lrm = at;
}
hrm = hrm + (t >> 16);
}
{
unsigned t = hxm * lym; /* 16 * 16 => 32 */
{
unsigned at = lrm + (t << 16);
hrm += at < lrm;
lrm = at;
}
hrm = hrm + (t >> 16);
}
rm = hrm | (lrm != 0);
}
/*
* sticky shift down to normal rounding precision
*/
if ((int) rm < 0) {
rm = (rm >> (32 - (SP_MBITS + 1 + 3))) |
((rm << (SP_MBITS + 1 + 3)) != 0);
re++;
} else {
rm = (rm >> (32 - (SP_MBITS + 1 + 3 + 1))) |
((rm << (SP_MBITS + 1 + 3 + 1)) != 0);
}
assert(rm & (SP_HIDDEN_BIT << 3));
SPNORMRET2(rs, re, rm, "mul", x, y);
}
}

58
arch/mips/math-emu/sp_scalb.c Normal file
View File

@@ -0,0 +1,58 @@
/* IEEE754 floating point arithmetic
* single precision
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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_scalb(ieee754sp x, int n)
{
COMPXSP;
CLEARCX;
EXPLODEXSP;
switch (xc) {
case IEEE754_CLASS_SNAN:
return ieee754sp_nanxcpt(x, "scalb", x, n);
case IEEE754_CLASS_QNAN:
case IEEE754_CLASS_INF:
case IEEE754_CLASS_ZERO:
return x;
case IEEE754_CLASS_DNORM:
SPDNORMX;
break;
case IEEE754_CLASS_NORM:
break;
}
SPNORMRET2(xs, xe + n, xm << 3, "scalb", x, n);
}
ieee754sp ieee754sp_ldexp(ieee754sp x, int n)
{
return ieee754sp_scalb(x, n);
}

84
arch/mips/math-emu/sp_simple.c Normal file
View File

@@ -0,0 +1,84 @@
/* IEEE754 floating point arithmetic
* single precision
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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_finite(ieee754sp x)
{
return SPBEXP(x) != SP_EMAX + 1 + SP_EBIAS;
}
ieee754sp ieee754sp_copysign(ieee754sp x, ieee754sp y)
{
CLEARCX;
SPSIGN(x) = SPSIGN(y);
return x;
}
ieee754sp ieee754sp_neg(ieee754sp x)
{
COMPXSP;
EXPLODEXSP;
CLEARCX;
FLUSHXSP;
if (xc == IEEE754_CLASS_SNAN) {
SETCX(IEEE754_INVALID_OPERATION);
return ieee754sp_nanxcpt(ieee754sp_indef(), "neg");
}
if (ieee754sp_isnan(x)) /* but not infinity */
return ieee754sp_nanxcpt(x, "neg", x);
/* quick fix up */
SPSIGN(x) ^= 1;
return x;
}
ieee754sp ieee754sp_abs(ieee754sp x)
{
COMPXSP;
EXPLODEXSP;
CLEARCX;
FLUSHXSP;
if (xc == IEEE754_CLASS_SNAN) {
SETCX(IEEE754_INVALID_OPERATION);
return ieee754sp_nanxcpt(ieee754sp_indef(), "abs");
}
if (ieee754sp_isnan(x)) /* but not infinity */
return ieee754sp_nanxcpt(x, "abs", x);
/* quick fix up */
SPSIGN(x) = 0;
return x;
}

117
arch/mips/math-emu/sp_sqrt.c Normal file
View File

@@ -0,0 +1,117 @@
/* IEEE754 floating point arithmetic
* single precision square root
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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_sqrt(ieee754sp x)
{
int ix, s, q, m, t, i;
unsigned int r;
COMPXSP;
/* take care of Inf and NaN */
EXPLODEXSP;
CLEARCX;
FLUSHXSP;
/* x == INF or NAN? */
switch (xc) {
case IEEE754_CLASS_QNAN:
/* sqrt(Nan) = Nan */
return ieee754sp_nanxcpt(x, "sqrt");
case IEEE754_CLASS_SNAN:
SETCX(IEEE754_INVALID_OPERATION);
return ieee754sp_nanxcpt(ieee754sp_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 ieee754sp_nanxcpt(ieee754sp_indef(), "sqrt");
}
/* sqrt(+Inf) = Inf */
return x;
case IEEE754_CLASS_DNORM:
case IEEE754_CLASS_NORM:
if (xs) {
/* sqrt(-x) = Nan */
SETCX(IEEE754_INVALID_OPERATION);
return ieee754sp_nanxcpt(ieee754sp_indef(), "sqrt");
}
break;
}
ix = x.bits;
/* normalize x */
m = (ix >> 23);
if (m == 0) { /* subnormal x */
for (i = 0; (ix & 0x00800000) == 0; i++)
ix <<= 1;
m -= i - 1;
}
m -= 127; /* unbias exponent */
ix = (ix & 0x007fffff) | 0x00800000;
if (m & 1) /* odd m, double x to make it even */
ix += ix;
m >>= 1; /* m = [m/2] */
/* generate sqrt(x) bit by bit */
ix += ix;
q = s = 0; /* q = sqrt(x) */
r = 0x01000000; /* r = moving bit from right to left */
while (r != 0) {
t = s + r;
if (t <= ix) {
s = t + r;
ix -= t;
q += r;
}
ix += ix;
r >>= 1;
}
if (ix != 0) {
SETCX(IEEE754_INEXACT);
switch (ieee754_csr.rm) {
case IEEE754_RP:
q += 2;
break;
case IEEE754_RN:
q += (q & 1);
break;
}
}
ix = (q >> 1) + 0x3f000000;
ix += (m << 23);
x.bits = ix;
return x;
}

184
arch/mips/math-emu/sp_sub.c Normal file
View File

@@ -0,0 +1,184 @@
/* IEEE754 floating point arithmetic
* single precision
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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_sub(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(), "sub", 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 handling
*/
case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
if (xs != ys)
return x;
SETCX(IEEE754_INVALID_OPERATION);
return ieee754sp_xcpt(ieee754sp_indef(), "sub", x, y);
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 ieee754sp_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 ieee754sp_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):
/* quick fix up */
DPSIGN(y) ^= 1;
return y;
case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
SPDNORMX;
case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
SPDNORMY;
break;
case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
SPDNORMX;
break;
case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
break;
}
/* flip sign of y and handle as add */
ys ^= 1;
assert(xm & SP_HIDDEN_BIT);
assert(ym & SP_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;
SPXSRSYn(s);
} else if (ye > xe) {
/* have to shift x fraction right to align
*/
int s = ye - xe;
SPXSRSXn(s);
}
assert(xe == ye);
assert(xe <= SP_EMAX);
if (xs == ys) {
/* generate 28 bit result of adding two 27 bit numbers
*/
xm = xm + ym;
xe = xe;
xs = xs;
if (xm >> (SP_MBITS + 1 + 3)) { /* carry out */
SPXSRSX1(); /* shift preserving sticky */
}
} else {
if (xm >= ym) {
xm = xm - ym;
xe = xe;
xs = xs;
} else {
xm = ym - xm;
xe = xe;
xs = ys;
}
if (xm == 0) {
if (ieee754_csr.rm == IEEE754_RD)
return ieee754sp_zero(1); /* round negative inf. => sign = -1 */
else
return ieee754sp_zero(0); /* other round modes => sign = 1 */
}
/* normalize to rounding precision
*/
while ((xm >> (SP_MBITS + 3)) == 0) {
xm <<= 1;
xe--;
}
}
SPNORMRET2(xs, xe, xm, "sub", x, y);
}

128
arch/mips/math-emu/sp_tint.c Normal file
View File

@@ -0,0 +1,128 @@
/* IEEE754 floating point arithmetic
* single precision
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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 <linux/kernel.h>
#include "ieee754sp.h"
int ieee754sp_tint(ieee754sp x)
{
COMPXSP;
CLEARCX;
EXPLODEXSP;
FLUSHXSP;
switch (xc) {
case IEEE754_CLASS_SNAN:
case IEEE754_CLASS_QNAN:
case IEEE754_CLASS_INF:
SETCX(IEEE754_INVALID_OPERATION);
return ieee754si_xcpt(ieee754si_indef(), "sp_tint", x);
case IEEE754_CLASS_ZERO:
return 0;
case IEEE754_CLASS_DNORM:
case IEEE754_CLASS_NORM:
break;
}
if (xe >= 31) {
/* look for valid corner case */
if (xe == 31 && xs && xm == SP_HIDDEN_BIT)
return -0x80000000;
/* Set invalid. We will only use overflow for floating
point overflow */
SETCX(IEEE754_INVALID_OPERATION);
return ieee754si_xcpt(ieee754si_indef(), "sp_tint", x);
}
/* oh gawd */
if (xe > SP_MBITS) {
xm <<= xe - SP_MBITS;
} else {
u32 residue;
int round;
int sticky;
int odd;
if (xe < -1) {
residue = xm;
round = 0;
sticky = residue != 0;
xm = 0;
}
else {
/* Shifting a u32 32 times does not work,
* so we do it in two steps. Be aware that xe
* may be -1 */
residue = xm << (xe + 1);
residue <<= 31 - SP_MBITS;
round = (residue >> 31) != 0;
sticky = (residue << 1) != 0;
xm >>= SP_MBITS - xe;
}
odd = (xm & 0x1) != 0x0;
switch (ieee754_csr.rm) {
case IEEE754_RN:
if (round && (sticky || odd))
xm++;
break;
case IEEE754_RZ:
break;
case IEEE754_RU: /* toward +Infinity */
if ((round || sticky) && !xs)
xm++;
break;
case IEEE754_RD: /* toward -Infinity */
if ((round || sticky) && xs)
xm++;
break;
}
if ((xm >> 31) != 0) {
/* This can happen after rounding */
SETCX(IEEE754_INVALID_OPERATION);
return ieee754si_xcpt(ieee754si_indef(), "sp_tint", x);
}
if (round || sticky)
SETCX(IEEE754_INEXACT);
}
if (xs)
return -xm;
else
return xm;
}
unsigned int ieee754sp_tuns(ieee754sp x)
{
ieee754sp hb = ieee754sp_1e31();
/* what if x < 0 ?? */
if (ieee754sp_lt(x, hb))
return (unsigned) ieee754sp_tint(x);
return (unsigned) ieee754sp_tint(ieee754sp_sub(x, hb)) |
((unsigned) 1 << 31);
}

123
arch/mips/math-emu/sp_tlong.c Normal file
View File

@@ -0,0 +1,123 @@
/* IEEE754 floating point arithmetic
* single precision
*/
/*
* MIPS floating point support
* Copyright (C) 1994-2000 Algorithmics Ltd.
* 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"
s64 ieee754sp_tlong(ieee754sp x)
{
COMPXDP; /* <-- need 64-bit mantissa tmp */
CLEARCX;
EXPLODEXSP;
FLUSHXSP;
switch (xc) {
case IEEE754_CLASS_SNAN:
case IEEE754_CLASS_QNAN:
case IEEE754_CLASS_INF:
SETCX(IEEE754_INVALID_OPERATION);
return ieee754di_xcpt(ieee754di_indef(), "sp_tlong", x);
case IEEE754_CLASS_ZERO:
return 0;
case IEEE754_CLASS_DNORM:
case IEEE754_CLASS_NORM:
break;
}
if (xe >= 63) {
/* look for valid corner case */
if (xe == 63 && xs && xm == SP_HIDDEN_BIT)
return -0x8000000000000000LL;
/* Set invalid. We will only use overflow for floating
point overflow */
SETCX(IEEE754_INVALID_OPERATION);
return ieee754di_xcpt(ieee754di_indef(), "sp_tlong", x);
}
/* oh gawd */
if (xe > SP_MBITS) {
xm <<= xe - SP_MBITS;
} else if (xe < SP_MBITS) {
u32 residue;
int round;
int sticky;
int odd;
if (xe < -1) {
residue = xm;
round = 0;
sticky = residue != 0;
xm = 0;
}
else {
residue = xm << (32 - SP_MBITS + xe);
round = (residue >> 31) != 0;
sticky = (residue << 1) != 0;
xm >>= SP_MBITS - xe;
}
odd = (xm & 0x1) != 0x0;
switch (ieee754_csr.rm) {
case IEEE754_RN:
if (round && (sticky || odd))
xm++;
break;
case IEEE754_RZ:
break;
case IEEE754_RU: /* toward +Infinity */
if ((round || sticky) && !xs)
xm++;
break;
case IEEE754_RD: /* toward -Infinity */
if ((round || sticky) && xs)
xm++;
break;
}
if ((xm >> 63) != 0) {
/* This can happen after rounding */
SETCX(IEEE754_INVALID_OPERATION);
return ieee754di_xcpt(ieee754di_indef(), "sp_tlong", x);
}
if (round || sticky)
SETCX(IEEE754_INEXACT);
}
if (xs)
return -xm;
else
return xm;
}
u64 ieee754sp_tulong(ieee754sp x)
{
ieee754sp hb = ieee754sp_1e63();
/* what if x < 0 ?? */
if (ieee754sp_lt(x, hb))
return (u64) ieee754sp_tlong(x);
return (u64) ieee754sp_tlong(ieee754sp_sub(x, hb)) |
(1ULL << 63);
}