2 * Linux/PA-RISC Project (http://www.parisc-linux.org/)
4 * Floating-point emulation code
5 * Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org>
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)
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.
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
25 * @(#) pa/spmath/sfadd.c $Revision: 1.1 $
28 * Single_add: add two single precision values.
30 * External Interfaces:
31 * sgl_fadd(leftptr, rightptr, dstptr, status)
33 * Internal Interfaces:
36 * <<please update with a overview of the operation of this file>>
43 #include "sgl_float.h"
46 * Single_add: add two single precision values.
50 sgl_floating_point
*leftptr
,
51 sgl_floating_point
*rightptr
,
52 sgl_floating_point
*dstptr
,
55 register unsigned int left
, right
, result
, extent
;
56 register unsigned int signless_upper_left
, signless_upper_right
, save
;
59 register int result_exponent
, right_exponent
, diff_exponent
;
60 register int sign_save
, jumpsize
;
61 register boolean inexact
= FALSE
;
62 register boolean underflowtrap
;
64 /* Create local copies of the numbers */
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
);
73 * check first operand for NaN's or infinity
75 if ((result_exponent
= Sgl_exponent(left
)) == SGL_INFINITY_EXPONENT
)
77 if (Sgl_iszero_mantissa(left
))
79 if (Sgl_isnotnan(right
))
81 if (Sgl_isinfinity(right
) && save
!=0)
84 * invalid since operands are opposite signed infinity's
86 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION
);
88 Sgl_makequietnan(result
);
102 * is NaN; signaling or quiet?
104 if (Sgl_isone_signaling(left
))
106 /* trap if INVALIDTRAP enabled */
107 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION
);
113 * is second operand a signaling NaN?
115 else if (Sgl_is_signalingnan(right
))
117 /* trap if INVALIDTRAP enabled */
118 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION
);
121 Sgl_set_quiet(right
);
131 } /* End left NaN or Infinity processing */
133 * check second operand for NaN's or infinity
135 if (Sgl_isinfinity_exponent(right
))
137 if (Sgl_iszero_mantissa(right
))
139 /* return infinity */
144 * is NaN; signaling or quiet?
146 if (Sgl_isone_signaling(right
))
148 /* trap if INVALIDTRAP enabled */
149 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION
);
152 Sgl_set_quiet(right
);
159 } /* End right NaN or Infinity processing */
161 /* Invariant: Must be dealing with finite numbers */
163 /* Compare operands by removing the sign */
164 Sgl_copytoint_exponentmantissa(left
,signless_upper_left
);
165 Sgl_copytoint_exponentmantissa(right
,signless_upper_right
);
167 /* sign difference selects add or sub operation. */
168 if(Sgl_ismagnitudeless(signless_upper_left
,signless_upper_right
))
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
);
176 /* Invariant: left is not smaller than right. */
178 if((right_exponent
= Sgl_exponent(right
)) == 0)
180 /* Denormalized operands. First look for zeroes */
181 if(Sgl_iszero_mantissa(right
))
184 if(Sgl_iszero_exponentmantissa(left
))
186 /* Both operands are zeros */
187 if(Is_rounding_mode(ROUNDMINUS
))
189 Sgl_or_signs(left
,/*with*/right
);
193 Sgl_and_signs(left
,/*with*/right
);
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() )
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
);
210 return(UNDERFLOWEXCEPTION
);
217 /* Neither are zeroes */
218 Sgl_clear_sign(right
); /* Exponent is already cleared */
219 if(result_exponent
== 0 )
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 )
226 Sgl_subtract(left
,/*minus*/right
,/*into*/result
);
227 if(Sgl_iszero_mantissa(result
))
229 if(Is_rounding_mode(ROUNDMINUS
))
231 Sgl_setone_sign(result
);
235 Sgl_setzero_sign(result
);
243 Sgl_addition(left
,right
,/*into*/result
);
244 if(Sgl_isone_hidden(result
))
250 if(Is_underflowtrap_enabled())
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
);
259 return(UNDERFLOWEXCEPTION
);
264 right_exponent
= 1; /* Set exponent to reflect different bias
265 * with denomalized numbers. */
269 Sgl_clear_signexponent_set_hidden(right
);
271 Sgl_clear_exponent_set_hidden(left
);
272 diff_exponent
= result_exponent
- right_exponent
;
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
280 if(diff_exponent
> SGL_THRESHOLD
)
282 diff_exponent
= SGL_THRESHOLD
;
285 /* Align right operand by shifting to right */
286 Sgl_right_align(/*operand*/right
,/*shifted by*/diff_exponent
,
287 /*and lower to*/extent
);
289 /* Treat sum and difference of the operands separately. */
290 if( (/*signed*/int) save
< 0 )
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.
297 Sgl_subtract_withextension(left
,/*minus*/right
,/*with*/extent
,/*into*/result
);
298 if(Sgl_iszero_hidden(result
))
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
);
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.
317 if(Sgl_iszero(result
))
318 /* Must have been "x-x" or "x+(-x)". */
320 if(Is_rounding_mode(ROUNDMINUS
)) Sgl_setone_sign(result
);
325 /* Look to see if normalization is finished. */
326 if(Sgl_isone_hidden(result
))
328 if(result_exponent
==0)
330 /* Denormalized, exponent should be zero. Left operand *
331 * was normalized, so extent (guard, round) was zero */
336 /* No further normalization is needed. */
337 Sgl_set_sign(result
,/*using*/sign_save
);
338 Ext_leftshiftby1(extent
);
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
;
348 /* Shift extension to complete one bit of normalization and
349 * update exponent. */
350 Ext_leftshiftby1(extent
);
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. */
357 while(Sgl_iszero_hiddenhigh7mantissa(result
))
359 Sgl_leftshiftby8(result
);
360 if((result_exponent
-= 8) <= 0 && !underflowtrap
)
363 /* Now narrow it down to the nibble */
364 if(Sgl_iszero_hiddenhigh3mantissa(result
))
366 /* The lower nibble contains the normalizing one */
367 Sgl_leftshiftby4(result
);
368 if((result_exponent
-= 4) <= 0 && !underflowtrap
)
371 /* Select case were first bit is set (already normalized)
372 * otherwise select the proper shift. */
373 if((jumpsize
= Sgl_hiddenhigh3mantissa(result
)) > 7)
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
);
382 Sgl_sethigh4bits(result
,/*using*/sign_save
);
387 Sgl_leftshiftby3(result
);
388 result_exponent
-= 3;
394 Sgl_leftshiftby2(result
);
395 result_exponent
-= 2;
403 Sgl_leftshiftby1(result
);
404 result_exponent
-= 1;
408 if(result_exponent
> 0)
410 Sgl_set_exponent(result
,/*using*/result_exponent
);
412 return(NOEXCEPTION
); /* Sign bit is already set */
414 /* Fixup potential underflows */
416 if(Is_underflowtrap_enabled())
418 Sgl_set_sign(result
,sign_save
);
419 Sgl_setwrapped_exponent(result
,result_exponent
,unfl
);
421 /* inexact = FALSE; */
422 return(UNDERFLOWEXCEPTION
);
425 * Since we cannot get an inexact denormalized result,
428 Sgl_right_align(result
,/*by*/(1-result_exponent
),extent
);
429 Sgl_clear_signexponent(result
);
430 Sgl_set_sign(result
,sign_save
);
433 } /* end if(hidden...)... */
434 /* Fall through and round */
435 } /* end if(save < 0)... */
439 Sgl_addition(left
,right
,/*to*/result
);
440 if(Sgl_isone_hiddenoverflow(result
))
442 /* Prenormalization required. */
443 Sgl_rightshiftby1_withextent(result
,extent
,extent
);
444 Sgl_arithrightshiftby1(result
);
446 } /* end if hiddenoverflow... */
447 } /* end else ...add magnitudes... */
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. */
454 if(Ext_isnotzero(extent
))
457 switch(Rounding_mode())
459 case ROUNDNEAREST
: /* The default. */
460 if(Ext_isone_sign(extent
))
462 /* at least 1/2 ulp */
463 if(Ext_isnotzero_lower(extent
) ||
464 Sgl_isone_lowmantissa(result
))
466 /* either exactly half way and odd or more than 1/2ulp */
467 Sgl_increment(result
);
473 if(Sgl_iszero_sign(result
))
475 /* Round up positive results */
476 Sgl_increment(result
);
481 if(Sgl_isone_sign(result
))
483 /* Round down negative results */
484 Sgl_increment(result
);
488 /* truncate is simple */
489 } /* end switch... */
490 if(Sgl_isone_hiddenoverflow(result
)) result_exponent
++;
492 if(result_exponent
== SGL_INFINITY_EXPONENT
)
495 if(Is_overflowtrap_enabled())
497 Sgl_setwrapped_exponent(result
,result_exponent
,ovfl
);
500 if (Is_inexacttrap_enabled())
501 return(OVERFLOWEXCEPTION
| INEXACTEXCEPTION
);
502 else Set_inexactflag();
503 return(OVERFLOWEXCEPTION
);
509 Sgl_setoverflow(result
);
512 else Sgl_set_exponent(result
,result_exponent
);
515 if(Is_inexacttrap_enabled()) return(INEXACTEXCEPTION
);
516 else Set_inexactflag();
This page took 0.040554 seconds and 5 git commands to generate.