Linux Kernel  3.7.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
sfadd.c
Go to the documentation of this file.
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) <[email protected]>
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
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  */
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 */
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 */
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 */
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 */
185  {
186  /* Both operands are zeros */
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  {
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  }
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  {
270  }
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  {
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 */
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:
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 */
496  {
497  Sgl_setwrapped_exponent(result,result_exponent,ovfl);
498  *dstptr = result;
499  if (inexact)
502  else Set_inexactflag();
503  return(OVERFLOWEXCEPTION);
504  }
505  else
506  {
508  inexact = TRUE;
509  Sgl_setoverflow(result);
510  }
511  }
512  else Sgl_set_exponent(result,result_exponent);
513  *dstptr = result;
514  if(inexact)
516  else Set_inexactflag();
517  return(NOEXCEPTION);
518  }