Linux-2.6.12-rc2
[linux-2.6.git] / arch / parisc / math-emu / dfdiv.c
1 /*
2  * Linux/PA-RISC Project (http://www.parisc-linux.org/)
3  *
4  * Floating-point emulation code
5  *  Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org>
6  *
7  *    This program is free software; you can redistribute it and/or modify
8  *    it under the terms of the GNU General Public License as published by
9  *    the Free Software Foundation; either version 2, or (at your option)
10  *    any later version.
11  *
12  *    This program is distributed in the hope that it will be useful,
13  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
14  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15  *    GNU General Public License for more details.
16  *
17  *    You should have received a copy of the GNU General Public License
18  *    along with this program; if not, write to the Free Software
19  *    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
20  */
21 /*
22  * BEGIN_DESC
23  *
24  *  File:
25  *      @(#)    pa/spmath/dfdiv.c               $Revision: 1.1 $
26  *
27  *  Purpose:
28  *      Double Precision Floating-point Divide
29  *
30  *  External Interfaces:
31  *      dbl_fdiv(srcptr1,srcptr2,dstptr,status)
32  *
33  *  Internal Interfaces:
34  *
35  *  Theory:
36  *      <<please update with a overview of the operation of this file>>
37  *
38  * END_DESC
39 */
40
41
42 #include "float.h"
43 #include "dbl_float.h"
44
45 /*
46  *  Double Precision Floating-point Divide
47  */
48
49 int
50 dbl_fdiv (dbl_floating_point * srcptr1, dbl_floating_point * srcptr2,
51           dbl_floating_point * dstptr, unsigned int *status)
52 {
53         register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
54         register unsigned int opnd3p1, opnd3p2, resultp1, resultp2;
55         register int dest_exponent, count;
56         register boolean inexact = FALSE, guardbit = FALSE, stickybit = FALSE;
57         boolean is_tiny;
58
59         Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
60         Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
61         /* 
62          * set sign bit of result 
63          */
64         if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1)) 
65                 Dbl_setnegativezerop1(resultp1);  
66         else Dbl_setzerop1(resultp1);
67         /*
68          * check first operand for NaN's or infinity
69          */
70         if (Dbl_isinfinity_exponent(opnd1p1)) {
71                 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
72                         if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
73                                 if (Dbl_isinfinity(opnd2p1,opnd2p2)) {
74                                         /* 
75                                          * invalid since both operands 
76                                          * are infinity 
77                                          */
78                                         if (Is_invalidtrap_enabled())
79                                                 return(INVALIDEXCEPTION);
80                                         Set_invalidflag();
81                                         Dbl_makequietnan(resultp1,resultp2);
82                                         Dbl_copytoptr(resultp1,resultp2,dstptr);
83                                         return(NOEXCEPTION);
84                                 }
85                                 /*
86                                  * return infinity
87                                  */
88                                 Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
89                                 Dbl_copytoptr(resultp1,resultp2,dstptr);
90                                 return(NOEXCEPTION);
91                         }
92                 }
93                 else {
94                         /*
95                          * is NaN; signaling or quiet?
96                          */
97                         if (Dbl_isone_signaling(opnd1p1)) {
98                                 /* trap if INVALIDTRAP enabled */
99                                 if (Is_invalidtrap_enabled())
100                                         return(INVALIDEXCEPTION);
101                                 /* make NaN quiet */
102                                 Set_invalidflag();
103                                 Dbl_set_quiet(opnd1p1);
104                         }
105                         /* 
106                          * is second operand a signaling NaN? 
107                          */
108                         else if (Dbl_is_signalingnan(opnd2p1)) {
109                                 /* trap if INVALIDTRAP enabled */
110                                 if (Is_invalidtrap_enabled())
111                                         return(INVALIDEXCEPTION);
112                                 /* make NaN quiet */
113                                 Set_invalidflag();
114                                 Dbl_set_quiet(opnd2p1);
115                                 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
116                                 return(NOEXCEPTION);
117                         }
118                         /*
119                          * return quiet NaN
120                          */
121                         Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
122                         return(NOEXCEPTION);
123                 }
124         }
125         /*
126          * check second operand for NaN's or infinity
127          */
128         if (Dbl_isinfinity_exponent(opnd2p1)) {
129                 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
130                         /*
131                          * return zero
132                          */
133                         Dbl_setzero_exponentmantissa(resultp1,resultp2);
134                         Dbl_copytoptr(resultp1,resultp2,dstptr);
135                         return(NOEXCEPTION);
136                 }
137                 /*
138                  * is NaN; signaling or quiet?
139                  */
140                 if (Dbl_isone_signaling(opnd2p1)) {
141                         /* trap if INVALIDTRAP enabled */
142                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
143                         /* make NaN quiet */
144                         Set_invalidflag();
145                         Dbl_set_quiet(opnd2p1);
146                 }
147                 /*
148                  * return quiet NaN
149                  */
150                 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
151                 return(NOEXCEPTION);
152         }
153         /*
154          * check for division by zero
155          */
156         if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
157                 if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) {
158                         /* invalid since both operands are zero */
159                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
160                         Set_invalidflag();
161                         Dbl_makequietnan(resultp1,resultp2);
162                         Dbl_copytoptr(resultp1,resultp2,dstptr);
163                         return(NOEXCEPTION);
164                 }
165                 if (Is_divisionbyzerotrap_enabled())
166                         return(DIVISIONBYZEROEXCEPTION);
167                 Set_divisionbyzeroflag();
168                 Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
169                 Dbl_copytoptr(resultp1,resultp2,dstptr);
170                 return(NOEXCEPTION);
171         }
172         /*
173          * Generate exponent 
174          */
175         dest_exponent = Dbl_exponent(opnd1p1) - Dbl_exponent(opnd2p1) + DBL_BIAS;
176
177         /*
178          * Generate mantissa
179          */
180         if (Dbl_isnotzero_exponent(opnd1p1)) {
181                 /* set hidden bit */
182                 Dbl_clear_signexponent_set_hidden(opnd1p1);
183         }
184         else {
185                 /* check for zero */
186                 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
187                         Dbl_setzero_exponentmantissa(resultp1,resultp2);
188                         Dbl_copytoptr(resultp1,resultp2,dstptr);
189                         return(NOEXCEPTION);
190                 }
191                 /* is denormalized, want to normalize */
192                 Dbl_clear_signexponent(opnd1p1);
193                 Dbl_leftshiftby1(opnd1p1,opnd1p2);
194                 Dbl_normalize(opnd1p1,opnd1p2,dest_exponent);
195         }
196         /* opnd2 needs to have hidden bit set with msb in hidden bit */
197         if (Dbl_isnotzero_exponent(opnd2p1)) {
198                 Dbl_clear_signexponent_set_hidden(opnd2p1);
199         }
200         else {
201                 /* is denormalized; want to normalize */
202                 Dbl_clear_signexponent(opnd2p1);
203                 Dbl_leftshiftby1(opnd2p1,opnd2p2);
204                 while (Dbl_iszero_hiddenhigh7mantissa(opnd2p1)) {
205                         dest_exponent+=8;
206                         Dbl_leftshiftby8(opnd2p1,opnd2p2);
207                 }
208                 if (Dbl_iszero_hiddenhigh3mantissa(opnd2p1)) {
209                         dest_exponent+=4;
210                         Dbl_leftshiftby4(opnd2p1,opnd2p2);
211                 }
212                 while (Dbl_iszero_hidden(opnd2p1)) {
213                         dest_exponent++;
214                         Dbl_leftshiftby1(opnd2p1,opnd2p2);
215                 }
216         }
217
218         /* Divide the source mantissas */
219
220         /* 
221          * A non-restoring divide algorithm is used.
222          */
223         Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
224         Dbl_setzero(opnd3p1,opnd3p2);
225         for (count=1; count <= DBL_P && (opnd1p1 || opnd1p2); count++) {
226                 Dbl_leftshiftby1(opnd1p1,opnd1p2);
227                 Dbl_leftshiftby1(opnd3p1,opnd3p2);
228                 if (Dbl_iszero_sign(opnd1p1)) {
229                         Dbl_setone_lowmantissap2(opnd3p2);
230                         Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
231                 }
232                 else {
233                         Twoword_add(opnd1p1, opnd1p2, opnd2p1, opnd2p2);
234                 }
235         }
236         if (count <= DBL_P) {
237                 Dbl_leftshiftby1(opnd3p1,opnd3p2);
238                 Dbl_setone_lowmantissap2(opnd3p2);
239                 Dbl_leftshift(opnd3p1,opnd3p2,(DBL_P-count));
240                 if (Dbl_iszero_hidden(opnd3p1)) {
241                         Dbl_leftshiftby1(opnd3p1,opnd3p2);
242                         dest_exponent--;
243                 }
244         }
245         else {
246                 if (Dbl_iszero_hidden(opnd3p1)) {
247                         /* need to get one more bit of result */
248                         Dbl_leftshiftby1(opnd1p1,opnd1p2);
249                         Dbl_leftshiftby1(opnd3p1,opnd3p2);
250                         if (Dbl_iszero_sign(opnd1p1)) {
251                                 Dbl_setone_lowmantissap2(opnd3p2);
252                                 Twoword_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
253                         }
254                         else {
255                                 Twoword_add(opnd1p1,opnd1p2,opnd2p1,opnd2p2);
256                         }
257                         dest_exponent--;
258                 }
259                 if (Dbl_iszero_sign(opnd1p1)) guardbit = TRUE;
260                 stickybit = Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2);
261         }
262         inexact = guardbit | stickybit;
263
264         /* 
265          * round result 
266          */
267         if (inexact && (dest_exponent > 0 || Is_underflowtrap_enabled())) {
268                 Dbl_clear_signexponent(opnd3p1);
269                 switch (Rounding_mode()) {
270                         case ROUNDPLUS: 
271                                 if (Dbl_iszero_sign(resultp1)) 
272                                         Dbl_increment(opnd3p1,opnd3p2);
273                                 break;
274                         case ROUNDMINUS: 
275                                 if (Dbl_isone_sign(resultp1)) 
276                                         Dbl_increment(opnd3p1,opnd3p2);
277                                 break;
278                         case ROUNDNEAREST:
279                                 if (guardbit && (stickybit || 
280                                     Dbl_isone_lowmantissap2(opnd3p2))) {
281                                         Dbl_increment(opnd3p1,opnd3p2);
282                                 }
283                 }
284                 if (Dbl_isone_hidden(opnd3p1)) dest_exponent++;
285         }
286         Dbl_set_mantissa(resultp1,resultp2,opnd3p1,opnd3p2);
287
288         /* 
289          * Test for overflow
290          */
291         if (dest_exponent >= DBL_INFINITY_EXPONENT) {
292                 /* trap if OVERFLOWTRAP enabled */
293                 if (Is_overflowtrap_enabled()) {
294                         /*
295                          * Adjust bias of result
296                          */
297                         Dbl_setwrapped_exponent(resultp1,dest_exponent,ovfl);
298                         Dbl_copytoptr(resultp1,resultp2,dstptr);
299                         if (inexact) 
300                             if (Is_inexacttrap_enabled())
301                                 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
302                             else Set_inexactflag();
303                         return(OVERFLOWEXCEPTION);
304                 }
305                 Set_overflowflag();
306                 /* set result to infinity or largest number */
307                 Dbl_setoverflow(resultp1,resultp2);
308                 inexact = TRUE;
309         }
310         /* 
311          * Test for underflow
312          */
313         else if (dest_exponent <= 0) {
314                 /* trap if UNDERFLOWTRAP enabled */
315                 if (Is_underflowtrap_enabled()) {
316                         /*
317                          * Adjust bias of result
318                          */
319                         Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
320                         Dbl_copytoptr(resultp1,resultp2,dstptr);
321                         if (inexact) 
322                             if (Is_inexacttrap_enabled())
323                                 return(UNDERFLOWEXCEPTION | INEXACTEXCEPTION);
324                             else Set_inexactflag();
325                         return(UNDERFLOWEXCEPTION);
326                 }
327
328                 /* Determine if should set underflow flag */
329                 is_tiny = TRUE;
330                 if (dest_exponent == 0 && inexact) {
331                         switch (Rounding_mode()) {
332                         case ROUNDPLUS: 
333                                 if (Dbl_iszero_sign(resultp1)) {
334                                         Dbl_increment(opnd3p1,opnd3p2);
335                                         if (Dbl_isone_hiddenoverflow(opnd3p1))
336                                             is_tiny = FALSE;
337                                         Dbl_decrement(opnd3p1,opnd3p2);
338                                 }
339                                 break;
340                         case ROUNDMINUS: 
341                                 if (Dbl_isone_sign(resultp1)) {
342                                         Dbl_increment(opnd3p1,opnd3p2);
343                                         if (Dbl_isone_hiddenoverflow(opnd3p1))
344                                             is_tiny = FALSE;
345                                         Dbl_decrement(opnd3p1,opnd3p2);
346                                 }
347                                 break;
348                         case ROUNDNEAREST:
349                                 if (guardbit && (stickybit || 
350                                     Dbl_isone_lowmantissap2(opnd3p2))) {
351                                         Dbl_increment(opnd3p1,opnd3p2);
352                                         if (Dbl_isone_hiddenoverflow(opnd3p1))
353                                             is_tiny = FALSE;
354                                         Dbl_decrement(opnd3p1,opnd3p2);
355                                 }
356                                 break;
357                         }
358                 }
359
360                 /*
361                  * denormalize result or set to signed zero
362                  */
363                 stickybit = inexact;
364                 Dbl_denormalize(opnd3p1,opnd3p2,dest_exponent,guardbit,
365                  stickybit,inexact);
366
367                 /* return rounded number */ 
368                 if (inexact) {
369                         switch (Rounding_mode()) {
370                         case ROUNDPLUS:
371                                 if (Dbl_iszero_sign(resultp1)) {
372                                         Dbl_increment(opnd3p1,opnd3p2);
373                                 }
374                                 break;
375                         case ROUNDMINUS: 
376                                 if (Dbl_isone_sign(resultp1)) {
377                                         Dbl_increment(opnd3p1,opnd3p2);
378                                 }
379                                 break;
380                         case ROUNDNEAREST:
381                                 if (guardbit && (stickybit || 
382                                     Dbl_isone_lowmantissap2(opnd3p2))) {
383                                         Dbl_increment(opnd3p1,opnd3p2);
384                                 }
385                                 break;
386                         }
387                         if (is_tiny) Set_underflowflag();
388                 }
389                 Dbl_set_exponentmantissa(resultp1,resultp2,opnd3p1,opnd3p2);
390         }
391         else Dbl_set_exponent(resultp1,dest_exponent);
392         Dbl_copytoptr(resultp1,resultp2,dstptr);
393
394         /* check for inexact */
395         if (inexact) {
396                 if (Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
397                 else Set_inexactflag();
398         }
399         return(NOEXCEPTION);
400 }