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