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