| 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_div(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(), "div", 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 | 		SETCX(IEEE754_INVALID_OPERATION); | 
 | 76 | 		return ieee754dp_xcpt(ieee754dp_indef(), "div", x, y); | 
 | 77 |  | 
 | 78 | 	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF): | 
 | 79 | 	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF): | 
 | 80 | 	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF): | 
 | 81 | 		return ieee754dp_zero(xs ^ ys); | 
 | 82 |  | 
 | 83 | 	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO): | 
 | 84 | 	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM): | 
 | 85 | 	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM): | 
 | 86 | 		return ieee754dp_inf(xs ^ ys); | 
 | 87 |  | 
 | 88 | 		/* Zero handling | 
 | 89 | 		 */ | 
 | 90 |  | 
 | 91 | 	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO): | 
 | 92 | 		SETCX(IEEE754_INVALID_OPERATION); | 
 | 93 | 		return ieee754dp_xcpt(ieee754dp_indef(), "div", x, y); | 
 | 94 |  | 
 | 95 | 	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO): | 
 | 96 | 	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO): | 
 | 97 | 		SETCX(IEEE754_ZERO_DIVIDE); | 
 | 98 | 		return ieee754dp_xcpt(ieee754dp_inf(xs ^ ys), "div", x, y); | 
 | 99 |  | 
 | 100 | 	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM): | 
 | 101 | 	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM): | 
 | 102 | 		return ieee754dp_zero(xs == ys ? 0 : 1); | 
 | 103 |  | 
 | 104 | 	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM): | 
 | 105 | 		DPDNORMX; | 
 | 106 |  | 
 | 107 | 	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM): | 
 | 108 | 		DPDNORMY; | 
 | 109 | 		break; | 
 | 110 |  | 
 | 111 | 	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM): | 
 | 112 | 		DPDNORMX; | 
 | 113 | 		break; | 
 | 114 |  | 
 | 115 | 	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM): | 
 | 116 | 		break; | 
 | 117 | 	} | 
 | 118 | 	assert(xm & DP_HIDDEN_BIT); | 
 | 119 | 	assert(ym & DP_HIDDEN_BIT); | 
 | 120 |  | 
 | 121 | 	/* provide rounding space */ | 
 | 122 | 	xm <<= 3; | 
 | 123 | 	ym <<= 3; | 
 | 124 |  | 
 | 125 | 	{ | 
 | 126 | 		/* now the dirty work */ | 
 | 127 |  | 
 | 128 | 		u64 rm = 0; | 
 | 129 | 		int re = xe - ye; | 
 | 130 | 		u64 bm; | 
 | 131 |  | 
 | 132 | 		for (bm = DP_MBIT(DP_MBITS + 2); bm; bm >>= 1) { | 
 | 133 | 			if (xm >= ym) { | 
 | 134 | 				xm -= ym; | 
 | 135 | 				rm |= bm; | 
 | 136 | 				if (xm == 0) | 
 | 137 | 					break; | 
 | 138 | 			} | 
 | 139 | 			xm <<= 1; | 
 | 140 | 		} | 
 | 141 | 		rm <<= 1; | 
 | 142 | 		if (xm) | 
 | 143 | 			rm |= 1;	/* have remainder, set sticky */ | 
 | 144 |  | 
 | 145 | 		assert(rm); | 
 | 146 |  | 
 | 147 | 		/* normalise rm to rounding precision ? | 
 | 148 | 		 */ | 
 | 149 | 		while ((rm >> (DP_MBITS + 3)) == 0) { | 
 | 150 | 			rm <<= 1; | 
 | 151 | 			re--; | 
 | 152 | 		} | 
 | 153 |  | 
 | 154 | 		DPNORMRET2(xs == ys ? 0 : 1, re, rm, "div", x, y); | 
 | 155 | 	} | 
 | 156 | } |