ARM: 6979/1: mach-vt8500: add forgotten irq_data conversion
[linux-2.6.git] / arch / parisc / math-emu / dfadd.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/dfadd.c               $Revision: 1.1 $
26  *
27  *  Purpose:
28  *      Double_add: add two double precision values.
29  *
30  *  External Interfaces:
31  *      dbl_fadd(leftptr, rightptr, 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_add: add two double precision values.
47  */
48 dbl_fadd(
49     dbl_floating_point *leftptr,
50     dbl_floating_point *rightptr,
51     dbl_floating_point *dstptr,
52     unsigned int *status)
53 {
54     register unsigned int signless_upper_left, signless_upper_right, save;
55     register unsigned int leftp1, leftp2, rightp1, rightp2, extent;
56     register unsigned int resultp1 = 0, resultp2 = 0;
57     
58     register int result_exponent, right_exponent, diff_exponent;
59     register int sign_save, jumpsize;
60     register boolean inexact = FALSE;
61     register boolean underflowtrap;
62         
63     /* Create local copies of the numbers */
64     Dbl_copyfromptr(leftptr,leftp1,leftp2);
65     Dbl_copyfromptr(rightptr,rightp1,rightp2);
66
67     /* A zero "save" helps discover equal operands (for later),  *
68      * and is used in swapping operands (if needed).             */
69     Dbl_xortointp1(leftp1,rightp1,/*to*/save);
70
71     /*
72      * check first operand for NaN's or infinity
73      */
74     if ((result_exponent = Dbl_exponent(leftp1)) == DBL_INFINITY_EXPONENT)
75         {
76         if (Dbl_iszero_mantissa(leftp1,leftp2)) 
77             {
78             if (Dbl_isnotnan(rightp1,rightp2)) 
79                 {
80                 if (Dbl_isinfinity(rightp1,rightp2) && save!=0) 
81                     {
82                     /* 
83                      * invalid since operands are opposite signed infinity's
84                      */
85                     if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
86                     Set_invalidflag();
87                     Dbl_makequietnan(resultp1,resultp2);
88                     Dbl_copytoptr(resultp1,resultp2,dstptr);
89                     return(NOEXCEPTION);
90                     }
91                 /*
92                  * return infinity
93                  */
94                 Dbl_copytoptr(leftp1,leftp2,dstptr);
95                 return(NOEXCEPTION);
96                 }
97             }
98         else 
99             {
100             /*
101              * is NaN; signaling or quiet?
102              */
103             if (Dbl_isone_signaling(leftp1)) 
104                 {
105                 /* trap if INVALIDTRAP enabled */
106                 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
107                 /* make NaN quiet */
108                 Set_invalidflag();
109                 Dbl_set_quiet(leftp1);
110                 }
111             /* 
112              * is second operand a signaling NaN? 
113              */
114             else if (Dbl_is_signalingnan(rightp1)) 
115                 {
116                 /* trap if INVALIDTRAP enabled */
117                 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
118                 /* make NaN quiet */
119                 Set_invalidflag();
120                 Dbl_set_quiet(rightp1);
121                 Dbl_copytoptr(rightp1,rightp2,dstptr);
122                 return(NOEXCEPTION);
123                 }
124             /*
125              * return quiet NaN
126              */
127             Dbl_copytoptr(leftp1,leftp2,dstptr);
128             return(NOEXCEPTION);
129             }
130         } /* End left NaN or Infinity processing */
131     /*
132      * check second operand for NaN's or infinity
133      */
134     if (Dbl_isinfinity_exponent(rightp1)) 
135         {
136         if (Dbl_iszero_mantissa(rightp1,rightp2)) 
137             {
138             /* return infinity */
139             Dbl_copytoptr(rightp1,rightp2,dstptr);
140             return(NOEXCEPTION);
141             }
142         /*
143          * is NaN; signaling or quiet?
144          */
145         if (Dbl_isone_signaling(rightp1)) 
146             {
147             /* trap if INVALIDTRAP enabled */
148             if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
149             /* make NaN quiet */
150             Set_invalidflag();
151             Dbl_set_quiet(rightp1);
152             }
153         /*
154          * return quiet NaN
155          */
156         Dbl_copytoptr(rightp1,rightp2,dstptr);
157         return(NOEXCEPTION);
158         } /* End right NaN or Infinity processing */
159
160     /* Invariant: Must be dealing with finite numbers */
161
162     /* Compare operands by removing the sign */
163     Dbl_copytoint_exponentmantissap1(leftp1,signless_upper_left);
164     Dbl_copytoint_exponentmantissap1(rightp1,signless_upper_right);
165
166     /* sign difference selects add or sub operation. */
167     if(Dbl_ismagnitudeless(leftp2,rightp2,signless_upper_left,signless_upper_right))
168         {
169         /* Set the left operand to the larger one by XOR swap *
170          *  First finish the first word using "save"          */
171         Dbl_xorfromintp1(save,rightp1,/*to*/rightp1);
172         Dbl_xorfromintp1(save,leftp1,/*to*/leftp1);
173         Dbl_swap_lower(leftp2,rightp2);
174         result_exponent = Dbl_exponent(leftp1);
175         }
176     /* Invariant:  left is not smaller than right. */ 
177
178     if((right_exponent = Dbl_exponent(rightp1)) == 0)
179         {
180         /* Denormalized operands.  First look for zeroes */
181         if(Dbl_iszero_mantissa(rightp1,rightp2)) 
182             {
183             /* right is zero */
184             if(Dbl_iszero_exponentmantissa(leftp1,leftp2))
185                 {
186                 /* Both operands are zeros */
187                 if(Is_rounding_mode(ROUNDMINUS))
188                     {
189                     Dbl_or_signs(leftp1,/*with*/rightp1);
190                     }
191                 else
192                     {
193                     Dbl_and_signs(leftp1,/*with*/rightp1);
194                     }
195                 }
196             else 
197                 {
198                 /* Left is not a zero and must be the result.  Trapped
199                  * underflows are signaled if left is denormalized.  Result
200                  * is always exact. */
201                 if( (result_exponent == 0) && Is_underflowtrap_enabled() )
202                     {
203                     /* need to normalize results mantissa */
204                     sign_save = Dbl_signextendedsign(leftp1);
205                     Dbl_leftshiftby1(leftp1,leftp2);
206                     Dbl_normalize(leftp1,leftp2,result_exponent);
207                     Dbl_set_sign(leftp1,/*using*/sign_save);
208                     Dbl_setwrapped_exponent(leftp1,result_exponent,unfl);
209                     Dbl_copytoptr(leftp1,leftp2,dstptr);
210                     /* inexact = FALSE */
211                     return(UNDERFLOWEXCEPTION);
212                     }
213                 }
214             Dbl_copytoptr(leftp1,leftp2,dstptr);
215             return(NOEXCEPTION);
216             }
217
218         /* Neither are zeroes */
219         Dbl_clear_sign(rightp1);        /* Exponent is already cleared */
220         if(result_exponent == 0 )
221             {
222             /* Both operands are denormalized.  The result must be exact
223              * and is simply calculated.  A sum could become normalized and a
224              * difference could cancel to a true zero. */
225             if( (/*signed*/int) save < 0 )
226                 {
227                 Dbl_subtract(leftp1,leftp2,/*minus*/rightp1,rightp2,
228                 /*into*/resultp1,resultp2);
229                 if(Dbl_iszero_mantissa(resultp1,resultp2))
230                     {
231                     if(Is_rounding_mode(ROUNDMINUS))
232                         {
233                         Dbl_setone_sign(resultp1);
234                         }
235                     else
236                         {
237                         Dbl_setzero_sign(resultp1);
238                         }
239                     Dbl_copytoptr(resultp1,resultp2,dstptr);
240                     return(NOEXCEPTION);
241                     }
242                 }
243             else
244                 {
245                 Dbl_addition(leftp1,leftp2,rightp1,rightp2,
246                 /*into*/resultp1,resultp2);
247                 if(Dbl_isone_hidden(resultp1))
248                     {
249                     Dbl_copytoptr(resultp1,resultp2,dstptr);
250                     return(NOEXCEPTION);
251                     }
252                 }
253             if(Is_underflowtrap_enabled())
254                 {
255                 /* need to normalize result */
256                 sign_save = Dbl_signextendedsign(resultp1);
257                 Dbl_leftshiftby1(resultp1,resultp2);
258                 Dbl_normalize(resultp1,resultp2,result_exponent);
259                 Dbl_set_sign(resultp1,/*using*/sign_save);
260                 Dbl_setwrapped_exponent(resultp1,result_exponent,unfl);
261                 Dbl_copytoptr(resultp1,resultp2,dstptr);
262                 /* inexact = FALSE */
263                 return(UNDERFLOWEXCEPTION);
264                 }
265             Dbl_copytoptr(resultp1,resultp2,dstptr);
266             return(NOEXCEPTION);
267             }
268         right_exponent = 1;     /* Set exponent to reflect different bias
269                                  * with denomalized numbers. */
270         }
271     else
272         {
273         Dbl_clear_signexponent_set_hidden(rightp1);
274         }
275     Dbl_clear_exponent_set_hidden(leftp1);
276     diff_exponent = result_exponent - right_exponent;
277
278     /* 
279      * Special case alignment of operands that would force alignment 
280      * beyond the extent of the extension.  A further optimization
281      * could special case this but only reduces the path length for this
282      * infrequent case.
283      */
284     if(diff_exponent > DBL_THRESHOLD)
285         {
286         diff_exponent = DBL_THRESHOLD;
287         }
288     
289     /* Align right operand by shifting to right */
290     Dbl_right_align(/*operand*/rightp1,rightp2,/*shifted by*/diff_exponent,
291     /*and lower to*/extent);
292
293     /* Treat sum and difference of the operands separately. */
294     if( (/*signed*/int) save < 0 )
295         {
296         /*
297          * Difference of the two operands.  Their can be no overflow.  A
298          * borrow can occur out of the hidden bit and force a post
299          * normalization phase.
300          */
301         Dbl_subtract_withextension(leftp1,leftp2,/*minus*/rightp1,rightp2,
302         /*with*/extent,/*into*/resultp1,resultp2);
303         if(Dbl_iszero_hidden(resultp1))
304             {
305             /* Handle normalization */
306             /* A straight forward algorithm would now shift the result
307              * and extension left until the hidden bit becomes one.  Not
308              * all of the extension bits need participate in the shift.
309              * Only the two most significant bits (round and guard) are
310              * needed.  If only a single shift is needed then the guard
311              * bit becomes a significant low order bit and the extension
312              * must participate in the rounding.  If more than a single 
313              * shift is needed, then all bits to the right of the guard 
314              * bit are zeros, and the guard bit may or may not be zero. */
315             sign_save = Dbl_signextendedsign(resultp1);
316             Dbl_leftshiftby1_withextent(resultp1,resultp2,extent,resultp1,resultp2);
317
318             /* Need to check for a zero result.  The sign and exponent
319              * fields have already been zeroed.  The more efficient test
320              * of the full object can be used.
321              */
322             if(Dbl_iszero(resultp1,resultp2))
323                 /* Must have been "x-x" or "x+(-x)". */
324                 {
325                 if(Is_rounding_mode(ROUNDMINUS)) Dbl_setone_sign(resultp1);
326                 Dbl_copytoptr(resultp1,resultp2,dstptr);
327                 return(NOEXCEPTION);
328                 }
329             result_exponent--;
330             /* Look to see if normalization is finished. */
331             if(Dbl_isone_hidden(resultp1))
332                 {
333                 if(result_exponent==0)
334                     {
335                     /* Denormalized, exponent should be zero.  Left operand *
336                      * was normalized, so extent (guard, round) was zero    */
337                     goto underflow;
338                     }
339                 else
340                     {
341                     /* No further normalization is needed. */
342                     Dbl_set_sign(resultp1,/*using*/sign_save);
343                     Ext_leftshiftby1(extent);
344                     goto round;
345                     }
346                 }
347
348             /* Check for denormalized, exponent should be zero.  Left    *
349              * operand was normalized, so extent (guard, round) was zero */
350             if(!(underflowtrap = Is_underflowtrap_enabled()) &&
351                result_exponent==0) goto underflow;
352
353             /* Shift extension to complete one bit of normalization and
354              * update exponent. */
355             Ext_leftshiftby1(extent);
356
357             /* Discover first one bit to determine shift amount.  Use a
358              * modified binary search.  We have already shifted the result
359              * one position right and still not found a one so the remainder
360              * of the extension must be zero and simplifies rounding. */
361             /* Scan bytes */
362             while(Dbl_iszero_hiddenhigh7mantissa(resultp1))
363                 {
364                 Dbl_leftshiftby8(resultp1,resultp2);
365                 if((result_exponent -= 8) <= 0  && !underflowtrap)
366                     goto underflow;
367                 }
368             /* Now narrow it down to the nibble */
369             if(Dbl_iszero_hiddenhigh3mantissa(resultp1))
370                 {
371                 /* The lower nibble contains the normalizing one */
372                 Dbl_leftshiftby4(resultp1,resultp2);
373                 if((result_exponent -= 4) <= 0 && !underflowtrap)
374                     goto underflow;
375                 }
376             /* Select case were first bit is set (already normalized)
377              * otherwise select the proper shift. */
378             if((jumpsize = Dbl_hiddenhigh3mantissa(resultp1)) > 7)
379                 {
380                 /* Already normalized */
381                 if(result_exponent <= 0) goto underflow;
382                 Dbl_set_sign(resultp1,/*using*/sign_save);
383                 Dbl_set_exponent(resultp1,/*using*/result_exponent);
384                 Dbl_copytoptr(resultp1,resultp2,dstptr);
385                 return(NOEXCEPTION);
386                 }
387             Dbl_sethigh4bits(resultp1,/*using*/sign_save);
388             switch(jumpsize) 
389                 {
390                 case 1:
391                     {
392                     Dbl_leftshiftby3(resultp1,resultp2);
393                     result_exponent -= 3;
394                     break;
395                     }
396                 case 2:
397                 case 3:
398                     {
399                     Dbl_leftshiftby2(resultp1,resultp2);
400                     result_exponent -= 2;
401                     break;
402                     }
403                 case 4:
404                 case 5:
405                 case 6:
406                 case 7:
407                     {
408                     Dbl_leftshiftby1(resultp1,resultp2);
409                     result_exponent -= 1;
410                     break;
411                     }
412                 }
413             if(result_exponent > 0) 
414                 {
415                 Dbl_set_exponent(resultp1,/*using*/result_exponent);
416                 Dbl_copytoptr(resultp1,resultp2,dstptr);
417                 return(NOEXCEPTION);    /* Sign bit is already set */
418                 }
419             /* Fixup potential underflows */
420           underflow:
421             if(Is_underflowtrap_enabled())
422                 {
423                 Dbl_set_sign(resultp1,sign_save);
424                 Dbl_setwrapped_exponent(resultp1,result_exponent,unfl);
425                 Dbl_copytoptr(resultp1,resultp2,dstptr);
426                 /* inexact = FALSE */
427                 return(UNDERFLOWEXCEPTION);
428                 }
429             /* 
430              * Since we cannot get an inexact denormalized result,
431              * we can now return.
432              */
433             Dbl_fix_overshift(resultp1,resultp2,(1-result_exponent),extent);
434             Dbl_clear_signexponent(resultp1);
435             Dbl_set_sign(resultp1,sign_save);
436             Dbl_copytoptr(resultp1,resultp2,dstptr);
437             return(NOEXCEPTION);
438             } /* end if(hidden...)... */
439         /* Fall through and round */
440         } /* end if(save < 0)... */
441     else 
442         {
443         /* Add magnitudes */
444         Dbl_addition(leftp1,leftp2,rightp1,rightp2,/*to*/resultp1,resultp2);
445         if(Dbl_isone_hiddenoverflow(resultp1))
446             {
447             /* Prenormalization required. */
448             Dbl_rightshiftby1_withextent(resultp2,extent,extent);
449             Dbl_arithrightshiftby1(resultp1,resultp2);
450             result_exponent++;
451             } /* end if hiddenoverflow... */
452         } /* end else ...add magnitudes... */
453     
454     /* Round the result.  If the extension is all zeros,then the result is
455      * exact.  Otherwise round in the correct direction.  No underflow is
456      * possible. If a postnormalization is necessary, then the mantissa is
457      * all zeros so no shift is needed. */
458   round:
459     if(Ext_isnotzero(extent))
460         {
461         inexact = TRUE;
462         switch(Rounding_mode())
463             {
464             case ROUNDNEAREST: /* The default. */
465             if(Ext_isone_sign(extent))
466                 {
467                 /* at least 1/2 ulp */
468                 if(Ext_isnotzero_lower(extent)  ||
469                   Dbl_isone_lowmantissap2(resultp2))
470                     {
471                     /* either exactly half way and odd or more than 1/2ulp */
472                     Dbl_increment(resultp1,resultp2);
473                     }
474                 }
475             break;
476
477             case ROUNDPLUS:
478             if(Dbl_iszero_sign(resultp1))
479                 {
480                 /* Round up positive results */
481                 Dbl_increment(resultp1,resultp2);
482                 }
483             break;
484             
485             case ROUNDMINUS:
486             if(Dbl_isone_sign(resultp1))
487                 {
488                 /* Round down negative results */
489                 Dbl_increment(resultp1,resultp2);
490                 }
491             
492             case ROUNDZERO:;
493             /* truncate is simple */
494             } /* end switch... */
495         if(Dbl_isone_hiddenoverflow(resultp1)) result_exponent++;
496         }
497     if(result_exponent == DBL_INFINITY_EXPONENT)
498         {
499         /* Overflow */
500         if(Is_overflowtrap_enabled())
501             {
502             Dbl_setwrapped_exponent(resultp1,result_exponent,ovfl);
503             Dbl_copytoptr(resultp1,resultp2,dstptr);
504             if (inexact)
505                 if (Is_inexacttrap_enabled())
506                         return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
507                 else Set_inexactflag();
508             return(OVERFLOWEXCEPTION);
509             }
510         else
511             {
512             inexact = TRUE;
513             Set_overflowflag();
514             Dbl_setoverflow(resultp1,resultp2);
515             }
516         }
517     else Dbl_set_exponent(resultp1,result_exponent);
518     Dbl_copytoptr(resultp1,resultp2,dstptr);
519     if(inexact) 
520         if(Is_inexacttrap_enabled())
521             return(INEXACTEXCEPTION);
522         else Set_inexactflag();
523     return(NOEXCEPTION);
524 }