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