1 /* IEEE754 floating point arithmetic
2  * single precision
3  */
4 /*
5  * MIPS floating point support
6  * Copyright (C) 1994-2000 Algorithmics Ltd.
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 "ieee754sp.h"
28 
ieee754sp_div(ieee754sp x,ieee754sp y)29 ieee754sp ieee754sp_div(ieee754sp x, ieee754sp y)
30 {
31 	COMPXSP;
32 	COMPYSP;
33 
34 	EXPLODEXSP;
35 	EXPLODEYSP;
36 
37 	CLEARCX;
38 
39 	FLUSHXSP;
40 	FLUSHYSP;
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 ieee754sp_nanxcpt(ieee754sp_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 ieee754sp_xcpt(ieee754sp_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 ieee754sp_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 ieee754sp_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 ieee754sp_xcpt(ieee754sp_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 ieee754sp_xcpt(ieee754sp_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 ieee754sp_zero(xs == ys ? 0 : 1);
103 
104 	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
105 		SPDNORMX;
106 
107 	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
108 		SPDNORMY;
109 		break;
110 
111 	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
112 		SPDNORMX;
113 		break;
114 
115 	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
116 		break;
117 	}
118 	assert(xm & SP_HIDDEN_BIT);
119 	assert(ym & SP_HIDDEN_BIT);
120 
121 	/* provide rounding space */
122 	xm <<= 3;
123 	ym <<= 3;
124 
125 	{
126 		/* now the dirty work */
127 
128 		unsigned rm = 0;
129 		int re = xe - ye;
130 		unsigned bm;
131 
132 		for (bm = SP_MBIT(SP_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 >> (SP_MBITS + 3)) == 0) {
150 			rm <<= 1;
151 			re--;
152 		}
153 
154 		SPNORMRET2(xs == ys ? 0 : 1, re, rm, "div", x, y);
155 	}
156 }
157