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 denormalized 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 | } |
506 | |