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