1 /* $NetBSD: dfadd.c,v 1.7 2023/07/31 20:48:04 andvar Exp $ */
2
3 /* $OpenBSD: dfadd.c,v 1.4 2001/03/29 03:58:17 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 * pmk1.1
27 */
28 /*
29 * (c) Copyright 1986 HEWLETT-PACKARD COMPANY
30 *
31 * To anyone who acknowledges that this file is provided "AS IS"
32 * without any express or implied warranty:
33 * permission to use, copy, modify, and distribute this file
34 * for any purpose is hereby granted without fee, provided that
35 * the above copyright notice and this notice appears in all
36 * copies, and that the name of Hewlett-Packard Company not be
37 * used in advertising or publicity pertaining to distribution
38 * of the software without specific, written prior permission.
39 * Hewlett-Packard Company makes no representations about the
40 * suitability of this software for any purpose.
41 */
42
43 #include <sys/cdefs.h>
44 __KERNEL_RCSID(0, "$NetBSD: dfadd.c,v 1.7 2023/07/31 20:48:04 andvar Exp $");
45
46 #include "../spmath/float.h"
47 #include "../spmath/dbl_float.h"
48
49 /*
50 * Double_add: add two double precision values.
51 */
52 int
dbl_fadd(dbl_floating_point * leftptr,dbl_floating_point * rightptr,dbl_floating_point * dstptr,unsigned int * status)53 dbl_fadd(dbl_floating_point *leftptr, dbl_floating_point *rightptr,
54 dbl_floating_point *dstptr, unsigned int *status)
55 {
56 register unsigned int signless_upper_left, signless_upper_right, save;
57 register unsigned int leftp1, leftp2, rightp1, rightp2, extent;
58 register unsigned int resultp1 = 0, resultp2 = 0;
59
60 register int result_exponent, right_exponent, diff_exponent;
61 register int sign_save, jumpsize;
62 register int inexact = false;
63 register int underflowtrap;
64
65 /* Create local copies of the numbers */
66 Dbl_copyfromptr(leftptr,leftp1,leftp2);
67 Dbl_copyfromptr(rightptr,rightp1,rightp2);
68
69 /* A zero "save" helps discover equal operands (for later), *
70 * and is used in swapping operands (if needed). */
71 Dbl_xortointp1(leftp1,rightp1,/*to*/save);
72
73 /*
74 * check first operand for NaN's or infinity
75 */
76 if ((result_exponent = Dbl_exponent(leftp1)) == DBL_INFINITY_EXPONENT)
77 {
78 if (Dbl_iszero_mantissa(leftp1,leftp2))
79 {
80 if (Dbl_isnotnan(rightp1,rightp2))
81 {
82 if (Dbl_isinfinity(rightp1,rightp2) && save!=0)
83 {
84 /*
85 * invalid since operands are opposite signed infinity's
86 */
87 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
88 Set_invalidflag();
89 Dbl_makequietnan(resultp1,resultp2);
90 Dbl_copytoptr(resultp1,resultp2,dstptr);
91 return(NOEXCEPTION);
92 }
93 /*
94 * return infinity
95 */
96 Dbl_copytoptr(leftp1,leftp2,dstptr);
97 return(NOEXCEPTION);
98 }
99 }
100 else
101 {
102 /*
103 * is NaN; signaling or quiet?
104 */
105 if (Dbl_isone_signaling(leftp1))
106 {
107 /* trap if INVALIDTRAP enabled */
108 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
109 /* make NaN quiet */
110 Set_invalidflag();
111 Dbl_set_quiet(leftp1);
112 }
113 /*
114 * is second operand a signaling NaN?
115 */
116 else if (Dbl_is_signalingnan(rightp1))
117 {
118 /* trap if INVALIDTRAP enabled */
119 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
120 /* make NaN quiet */
121 Set_invalidflag();
122 Dbl_set_quiet(rightp1);
123 Dbl_copytoptr(rightp1,rightp2,dstptr);
124 return(NOEXCEPTION);
125 }
126 /*
127 * return quiet NaN
128 */
129 Dbl_copytoptr(leftp1,leftp2,dstptr);
130 return(NOEXCEPTION);
131 }
132 } /* End left NaN or Infinity processing */
133 /*
134 * check second operand for NaN's or infinity
135 */
136 if (Dbl_isinfinity_exponent(rightp1))
137 {
138 if (Dbl_iszero_mantissa(rightp1,rightp2))
139 {
140 /* return infinity */
141 Dbl_copytoptr(rightp1,rightp2,dstptr);
142 return(NOEXCEPTION);
143 }
144 /*
145 * is NaN; signaling or quiet?
146 */
147 if (Dbl_isone_signaling(rightp1))
148 {
149 /* trap if INVALIDTRAP enabled */
150 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
151 /* make NaN quiet */
152 Set_invalidflag();
153 Dbl_set_quiet(rightp1);
154 }
155 /*
156 * return quiet NaN
157 */
158 Dbl_copytoptr(rightp1,rightp2,dstptr);
159 return(NOEXCEPTION);
160 } /* End right NaN or Infinity processing */
161
162 /* Invariant: Must be dealing with finite numbers */
163
164 /* Compare operands by removing the sign */
165 Dbl_copytoint_exponentmantissap1(leftp1,signless_upper_left);
166 Dbl_copytoint_exponentmantissap1(rightp1,signless_upper_right);
167
168 /* sign difference selects add or sub operation. */
169 if(Dbl_ismagnitudeless(leftp2,rightp2,signless_upper_left,signless_upper_right))
170 {
171 /* Set the left operand to the larger one by XOR swap *
172 * First finish the first word using "save" */
173 Dbl_xorfromintp1(save,rightp1,/*to*/rightp1);
174 Dbl_xorfromintp1(save,leftp1,/*to*/leftp1);
175 Dbl_swap_lower(leftp2,rightp2);
176 result_exponent = Dbl_exponent(leftp1);
177 }
178 /* Invariant: left is not smaller than right. */
179
180 if((right_exponent = Dbl_exponent(rightp1)) == 0)
181 {
182 /* Denormalized operands. First look for zeroes */
183 if(Dbl_iszero_mantissa(rightp1,rightp2))
184 {
185 /* right is zero */
186 if(Dbl_iszero_exponentmantissa(leftp1,leftp2))
187 {
188 /* Both operands are zeros */
189 if(Is_rounding_mode(ROUNDMINUS))
190 {
191 Dbl_or_signs(leftp1,/*with*/rightp1);
192 }
193 else
194 {
195 Dbl_and_signs(leftp1,/*with*/rightp1);
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 = Dbl_signextendedsign(leftp1);
207 Dbl_leftshiftby1(leftp1,leftp2);
208 Dbl_normalize(leftp1,leftp2,result_exponent);
209 Dbl_set_sign(leftp1,/*using*/sign_save);
210 Dbl_setwrapped_exponent(leftp1,result_exponent,unfl);
211 Dbl_copytoptr(leftp1,leftp2,dstptr);
212 /* inexact = false */
213 return(UNDERFLOWEXCEPTION);
214 }
215 }
216 Dbl_copytoptr(leftp1,leftp2,dstptr);
217 return(NOEXCEPTION);
218 }
219
220 /* Neither are zeroes */
221 Dbl_clear_sign(rightp1); /* Exponent is already cleared */
222 if(result_exponent == 0 )
223 {
224 /* Both operands are denormalized. The result must be exact
225 * and is simply calculated. A sum could become normalized and a
226 * difference could cancel to a true zero. */
227 if( (/*signed*/int) save < 0 )
228 {
229 Dbl_subtract(leftp1,leftp2,/*minus*/rightp1,rightp2,
230 /*into*/resultp1,resultp2);
231 if(Dbl_iszero_mantissa(resultp1,resultp2))
232 {
233 if(Is_rounding_mode(ROUNDMINUS))
234 {
235 Dbl_setone_sign(resultp1);
236 }
237 else
238 {
239 Dbl_setzero_sign(resultp1);
240 }
241 Dbl_copytoptr(resultp1,resultp2,dstptr);
242 return(NOEXCEPTION);
243 }
244 }
245 else
246 {
247 Dbl_addition(leftp1,leftp2,rightp1,rightp2,
248 /*into*/resultp1,resultp2);
249 if(Dbl_isone_hidden(resultp1))
250 {
251 Dbl_copytoptr(resultp1,resultp2,dstptr);
252 return(NOEXCEPTION);
253 }
254 }
255 if(Is_underflowtrap_enabled())
256 {
257 /* need to normalize result */
258 sign_save = Dbl_signextendedsign(resultp1);
259 Dbl_leftshiftby1(resultp1,resultp2);
260 Dbl_normalize(resultp1,resultp2,result_exponent);
261 Dbl_set_sign(resultp1,/*using*/sign_save);
262 Dbl_setwrapped_exponent(resultp1,result_exponent,unfl);
263 Dbl_copytoptr(resultp1,resultp2,dstptr);
264 /* inexact = false */
265 return(UNDERFLOWEXCEPTION);
266 }
267 Dbl_copytoptr(resultp1,resultp2,dstptr);
268 return(NOEXCEPTION);
269 }
270 right_exponent = 1; /* Set exponent to reflect different bias
271 * with denormalized numbers. */
272 }
273 else
274 {
275 Dbl_clear_signexponent_set_hidden(rightp1);
276 }
277 Dbl_clear_exponent_set_hidden(leftp1);
278 diff_exponent = result_exponent - right_exponent;
279
280 /*
281 * Special case alignment of operands that would force alignment
282 * beyond the extent of the extension. A further optimization
283 * could special case this but only reduces the path length for this
284 * infrequent case.
285 */
286 if(diff_exponent > DBL_THRESHOLD)
287 {
288 diff_exponent = DBL_THRESHOLD;
289 }
290
291 /* Align right operand by shifting to right */
292 Dbl_right_align(/*operand*/rightp1,rightp2,/*shifted by*/diff_exponent,
293 /*and lower to*/extent);
294
295 /* Treat sum and difference of the operands separately. */
296 if( (/*signed*/int) save < 0 )
297 {
298 /*
299 * Difference of the two operands. Their can be no overflow. A
300 * borrow can occur out of the hidden bit and force a post
301 * normalization phase.
302 */
303 Dbl_subtract_withextension(leftp1,leftp2,/*minus*/rightp1,rightp2,
304 /*with*/extent,/*into*/resultp1,resultp2);
305 if(Dbl_iszero_hidden(resultp1))
306 {
307 /* Handle normalization */
308 /* A straight forward algorithm would now shift the result
309 * and extension left until the hidden bit becomes one. Not
310 * all of the extension bits need participate in the shift.
311 * Only the two most significant bits (round and guard) are
312 * needed. If only a single shift is needed then the guard
313 * bit becomes a significant low order bit and the extension
314 * must participate in the rounding. If more than a single
315 * shift is needed, then all bits to the right of the guard
316 * bit are zeros, and the guard bit may or may not be zero. */
317 sign_save = Dbl_signextendedsign(resultp1);
318 Dbl_leftshiftby1_withextent(resultp1,resultp2,extent,resultp1,resultp2);
319
320 /* Need to check for a zero result. The sign and exponent
321 * fields have already been zeroed. The more efficient test
322 * of the full object can be used.
323 */
324 if(Dbl_iszero(resultp1,resultp2))
325 /* Must have been "x-x" or "x+(-x)". */
326 {
327 if(Is_rounding_mode(ROUNDMINUS)) Dbl_setone_sign(resultp1);
328 Dbl_copytoptr(resultp1,resultp2,dstptr);
329 return(NOEXCEPTION);
330 }
331 result_exponent--;
332 /* Look to see if normalization is finished. */
333 if(Dbl_isone_hidden(resultp1))
334 {
335 if(result_exponent==0)
336 {
337 /* Denormalized, exponent should be zero. Left operand *
338 * was normalized, so extent (guard, round) was zero */
339 goto underflow;
340 }
341 else
342 {
343 /* No further normalization is needed. */
344 Dbl_set_sign(resultp1,/*using*/sign_save);
345 Ext_leftshiftby1(extent);
346 goto round;
347 }
348 }
349
350 /* Check for denormalized, exponent should be zero. Left *
351 * operand was normalized, so extent (guard, round) was zero */
352 if(!(underflowtrap = Is_underflowtrap_enabled()) &&
353 result_exponent==0) goto underflow;
354
355 /* Shift extension to complete one bit of normalization and
356 * update exponent. */
357 Ext_leftshiftby1(extent);
358
359 /* Discover first one bit to determine shift amount. Use a
360 * modified binary search. We have already shifted the result
361 * one position right and still not found a one so the remainder
362 * of the extension must be zero and simplifies rounding. */
363 /* Scan bytes */
364 while(Dbl_iszero_hiddenhigh7mantissa(resultp1))
365 {
366 Dbl_leftshiftby8(resultp1,resultp2);
367 if((result_exponent -= 8) <= 0 && !underflowtrap)
368 goto underflow;
369 }
370 /* Now narrow it down to the nibble */
371 if(Dbl_iszero_hiddenhigh3mantissa(resultp1))
372 {
373 /* The lower nibble contains the normalizing one */
374 Dbl_leftshiftby4(resultp1,resultp2);
375 if((result_exponent -= 4) <= 0 && !underflowtrap)
376 goto underflow;
377 }
378 /* Select case were first bit is set (already normalized)
379 * otherwise select the proper shift. */
380 if((jumpsize = Dbl_hiddenhigh3mantissa(resultp1)) > 7)
381 {
382 /* Already normalized */
383 if(result_exponent <= 0) goto underflow;
384 Dbl_set_sign(resultp1,/*using*/sign_save);
385 Dbl_set_exponent(resultp1,/*using*/result_exponent);
386 Dbl_copytoptr(resultp1,resultp2,dstptr);
387 return(NOEXCEPTION);
388 }
389 Dbl_sethigh4bits(resultp1,/*using*/sign_save);
390 switch(jumpsize)
391 {
392 case 1:
393 {
394 Dbl_leftshiftby3(resultp1,resultp2);
395 result_exponent -= 3;
396 break;
397 }
398 case 2:
399 case 3:
400 {
401 Dbl_leftshiftby2(resultp1,resultp2);
402 result_exponent -= 2;
403 break;
404 }
405 case 4:
406 case 5:
407 case 6:
408 case 7:
409 {
410 Dbl_leftshiftby1(resultp1,resultp2);
411 result_exponent -= 1;
412 break;
413 }
414 }
415 if(result_exponent > 0)
416 {
417 Dbl_set_exponent(resultp1,/*using*/result_exponent);
418 Dbl_copytoptr(resultp1,resultp2,dstptr);
419 return(NOEXCEPTION); /* Sign bit is already set */
420 }
421 /* Fixup potential underflows */
422 underflow:
423 if(Is_underflowtrap_enabled())
424 {
425 Dbl_set_sign(resultp1,sign_save);
426 Dbl_setwrapped_exponent(resultp1,result_exponent,unfl);
427 Dbl_copytoptr(resultp1,resultp2,dstptr);
428 /* inexact = false */
429 return(UNDERFLOWEXCEPTION);
430 }
431 /*
432 * Since we cannot get an inexact denormalized result,
433 * we can now return.
434 */
435 Dbl_fix_overshift(resultp1,resultp2,(1-result_exponent),extent);
436 Dbl_clear_signexponent(resultp1);
437 Dbl_set_sign(resultp1,sign_save);
438 Dbl_copytoptr(resultp1,resultp2,dstptr);
439 return(NOEXCEPTION);
440 } /* end if(hidden...)... */
441 /* Fall through and round */
442 } /* end if(save < 0)... */
443 else
444 {
445 /* Add magnitudes */
446 Dbl_addition(leftp1,leftp2,rightp1,rightp2,/*to*/resultp1,resultp2);
447 if(Dbl_isone_hiddenoverflow(resultp1))
448 {
449 /* Prenormalization required. */
450 Dbl_rightshiftby1_withextent(resultp2,extent,extent);
451 Dbl_arithrightshiftby1(resultp1,resultp2);
452 result_exponent++;
453 } /* end if hiddenoverflow... */
454 } /* end else ...add magnitudes... */
455
456 /* Round the result. If the extension is all zeros,then the result is
457 * exact. Otherwise round in the correct direction. No underflow is
458 * possible. If a postnormalization is necessary, then the mantissa is
459 * all zeros so no shift is needed. */
460 round:
461 if(Ext_isnotzero(extent))
462 {
463 inexact = true;
464 switch(Rounding_mode())
465 {
466 case ROUNDNEAREST: /* The default. */
467 if(Ext_isone_sign(extent))
468 {
469 /* at least 1/2 ulp */
470 if(Ext_isnotzero_lower(extent) ||
471 Dbl_isone_lowmantissap2(resultp2))
472 {
473 /* either exactly half way and odd or more than 1/2ulp */
474 Dbl_increment(resultp1,resultp2);
475 }
476 }
477 break;
478
479 case ROUNDPLUS:
480 if(Dbl_iszero_sign(resultp1))
481 {
482 /* Round up positive results */
483 Dbl_increment(resultp1,resultp2);
484 }
485 break;
486
487 case ROUNDMINUS:
488 if(Dbl_isone_sign(resultp1))
489 {
490 /* Round down negative results */
491 Dbl_increment(resultp1,resultp2);
492 }
493
494 case ROUNDZERO:;
495 /* truncate is simple */
496 } /* end switch... */
497 if(Dbl_isone_hiddenoverflow(resultp1)) result_exponent++;
498 }
499 if(result_exponent == DBL_INFINITY_EXPONENT)
500 {
501 /* Overflow */
502 if(Is_overflowtrap_enabled())
503 {
504 Dbl_setwrapped_exponent(resultp1,result_exponent,ovfl);
505 Dbl_copytoptr(resultp1,resultp2,dstptr);
506 if (inexact) {
507 if (Is_inexacttrap_enabled())
508 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
509 else
510 Set_inexactflag();
511 }
512 return(OVERFLOWEXCEPTION);
513 }
514 else
515 {
516 inexact = true;
517 Set_overflowflag();
518 Dbl_setoverflow(resultp1,resultp2);
519 }
520 }
521 else Dbl_set_exponent(resultp1,result_exponent);
522 Dbl_copytoptr(resultp1,resultp2,dstptr);
523 if(inexact) {
524 if(Is_inexacttrap_enabled())
525 return(INEXACTEXCEPTION);
526 else
527 Set_inexactflag();
528 }
529 return(NOEXCEPTION);
530 }
531