1 /* $NetBSD: dfsub.c,v 1.7 2023/07/31 20:48:04 andvar Exp $ */
2
3 /* $OpenBSD: dfsub.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 /*
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: dfsub.c,v 1.7 2023/07/31 20:48:04 andvar Exp $");
46
47 #include "../spmath/float.h"
48 #include "../spmath/dbl_float.h"
49
50 /*
51 * Double_subtract: subtract two double precision values.
52 */
53 int
dbl_fsub(dbl_floating_point * leftptr,dbl_floating_point * rightptr,dbl_floating_point * dstptr,unsigned int * status)54 dbl_fsub(dbl_floating_point *leftptr, dbl_floating_point *rightptr,
55 dbl_floating_point *dstptr, unsigned int *status)
56 {
57 register unsigned int signless_upper_left, signless_upper_right, save;
58 register unsigned int leftp1, leftp2, rightp1, rightp2, extent;
59 register unsigned int resultp1 = 0, resultp2 = 0;
60
61 register int result_exponent, right_exponent, diff_exponent;
62 register int sign_save, jumpsize;
63 register int inexact = false, 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 same 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_invert_sign(rightp1);
142 Dbl_copytoptr(rightp1,rightp2,dstptr);
143 return(NOEXCEPTION);
144 }
145 /*
146 * is NaN; signaling or quiet?
147 */
148 if (Dbl_isone_signaling(rightp1))
149 {
150 /* trap if INVALIDTRAP enabled */
151 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
152 /* make NaN quiet */
153 Set_invalidflag();
154 Dbl_set_quiet(rightp1);
155 }
156 /*
157 * return quiet NaN
158 */
159 Dbl_copytoptr(rightp1,rightp2,dstptr);
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 Dbl_copytoint_exponentmantissap1(leftp1,signless_upper_left);
167 Dbl_copytoint_exponentmantissap1(rightp1,signless_upper_right);
168
169 /* sign difference selects add or sub operation. */
170 if(Dbl_ismagnitudeless(leftp2,rightp2,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 Dbl_xorfromintp1(save,rightp1,/*to*/rightp1);
175 Dbl_xorfromintp1(save,leftp1,/*to*/leftp1);
176 Dbl_swap_lower(leftp2,rightp2);
177 result_exponent = Dbl_exponent(leftp1);
178 Dbl_invert_sign(leftp1);
179 }
180 /* Invariant: left is not smaller than right. */
181
182 if((right_exponent = Dbl_exponent(rightp1)) == 0)
183 {
184 /* Denormalized operands. First look for zeroes */
185 if(Dbl_iszero_mantissa(rightp1,rightp2))
186 {
187 /* right is zero */
188 if(Dbl_iszero_exponentmantissa(leftp1,leftp2))
189 {
190 /* Both operands are zeros */
191 Dbl_invert_sign(rightp1);
192 if(Is_rounding_mode(ROUNDMINUS))
193 {
194 Dbl_or_signs(leftp1,/*with*/rightp1);
195 }
196 else
197 {
198 Dbl_and_signs(leftp1,/*with*/rightp1);
199 }
200 }
201 else
202 {
203 /* Left is not a zero and must be the result. Trapped
204 * underflows are signaled if left is denormalized. Result
205 * is always exact. */
206 if( (result_exponent == 0) && Is_underflowtrap_enabled() )
207 {
208 /* need to normalize results mantissa */
209 sign_save = Dbl_signextendedsign(leftp1);
210 Dbl_leftshiftby1(leftp1,leftp2);
211 Dbl_normalize(leftp1,leftp2,result_exponent);
212 Dbl_set_sign(leftp1,/*using*/sign_save);
213 Dbl_setwrapped_exponent(leftp1,result_exponent,unfl);
214 Dbl_copytoptr(leftp1,leftp2,dstptr);
215 /* inexact = false */
216 return(UNDERFLOWEXCEPTION);
217 }
218 }
219 Dbl_copytoptr(leftp1,leftp2,dstptr);
220 return(NOEXCEPTION);
221 }
222
223 /* Neither are zeroes */
224 Dbl_clear_sign(rightp1); /* Exponent is already cleared */
225 if(result_exponent == 0 )
226 {
227 /* Both operands are denormalized. The result must be exact
228 * and is simply calculated. A sum could become normalized and a
229 * difference could cancel to a true zero. */
230 if( (/*signed*/int) save >= 0 )
231 {
232 Dbl_subtract(leftp1,leftp2,/*minus*/rightp1,rightp2,
233 /*into*/resultp1,resultp2);
234 if(Dbl_iszero_mantissa(resultp1,resultp2))
235 {
236 if(Is_rounding_mode(ROUNDMINUS))
237 {
238 Dbl_setone_sign(resultp1);
239 }
240 else
241 {
242 Dbl_setzero_sign(resultp1);
243 }
244 Dbl_copytoptr(resultp1,resultp2,dstptr);
245 return(NOEXCEPTION);
246 }
247 }
248 else
249 {
250 Dbl_addition(leftp1,leftp2,rightp1,rightp2,
251 /*into*/resultp1,resultp2);
252 if(Dbl_isone_hidden(resultp1))
253 {
254 Dbl_copytoptr(resultp1,resultp2,dstptr);
255 return(NOEXCEPTION);
256 }
257 }
258 if(Is_underflowtrap_enabled())
259 {
260 /* need to normalize result */
261 sign_save = Dbl_signextendedsign(resultp1);
262 Dbl_leftshiftby1(resultp1,resultp2);
263 Dbl_normalize(resultp1,resultp2,result_exponent);
264 Dbl_set_sign(resultp1,/*using*/sign_save);
265 Dbl_setwrapped_exponent(resultp1,result_exponent,unfl);
266 Dbl_copytoptr(resultp1,resultp2,dstptr);
267 /* inexact = false */
268 return(UNDERFLOWEXCEPTION);
269 }
270 Dbl_copytoptr(resultp1,resultp2,dstptr);
271 return(NOEXCEPTION);
272 }
273 right_exponent = 1; /* Set exponent to reflect different bias
274 * with denormalized numbers. */
275 }
276 else
277 {
278 Dbl_clear_signexponent_set_hidden(rightp1);
279 }
280 Dbl_clear_exponent_set_hidden(leftp1);
281 diff_exponent = result_exponent - right_exponent;
282
283 /*
284 * Special case alignment of operands that would force alignment
285 * beyond the extent of the extension. A further optimization
286 * could special case this but only reduces the path length for this
287 * infrequent case.
288 */
289 if(diff_exponent > DBL_THRESHOLD)
290 {
291 diff_exponent = DBL_THRESHOLD;
292 }
293
294 /* Align right operand by shifting to right */
295 Dbl_right_align(/*operand*/rightp1,rightp2,/*shifted by*/diff_exponent,
296 /*and lower to*/extent);
297
298 /* Treat sum and difference of the operands separately. */
299 if( (/*signed*/int) save >= 0 )
300 {
301 /*
302 * Difference of the two operands. Their can be no overflow. A
303 * borrow can occur out of the hidden bit and force a post
304 * normalization phase.
305 */
306 Dbl_subtract_withextension(leftp1,leftp2,/*minus*/rightp1,rightp2,
307 /*with*/extent,/*into*/resultp1,resultp2);
308 if(Dbl_iszero_hidden(resultp1))
309 {
310 /* Handle normalization */
311 /* A straight forward algorithm would now shift the result
312 * and extension left until the hidden bit becomes one. Not
313 * all of the extension bits need participate in the shift.
314 * Only the two most significant bits (round and guard) are
315 * needed. If only a single shift is needed then the guard
316 * bit becomes a significant low order bit and the extension
317 * must participate in the rounding. If more than a single
318 * shift is needed, then all bits to the right of the guard
319 * bit are zeros, and the guard bit may or may not be zero. */
320 sign_save = Dbl_signextendedsign(resultp1);
321 Dbl_leftshiftby1_withextent(resultp1,resultp2,extent,resultp1,resultp2);
322
323 /* Need to check for a zero result. The sign and exponent
324 * fields have already been zeroed. The more efficient test
325 * of the full object can be used.
326 */
327 if(Dbl_iszero(resultp1,resultp2))
328 /* Must have been "x-x" or "x+(-x)". */
329 {
330 if(Is_rounding_mode(ROUNDMINUS)) Dbl_setone_sign(resultp1);
331 Dbl_copytoptr(resultp1,resultp2,dstptr);
332 return(NOEXCEPTION);
333 }
334 result_exponent--;
335 /* Look to see if normalization is finished. */
336 if(Dbl_isone_hidden(resultp1)) {
337 if(result_exponent==0) {
338 /* Denormalized, exponent should be zero. Left operand *
339 * was normalized, so extent (guard, round) was zero */
340 goto underflow;
341 } else {
342 /* No further normalization is needed. */
343 Dbl_set_sign(resultp1,/*using*/sign_save);
344 Ext_leftshiftby1(extent);
345 goto round;
346 }
347 }
348
349 /* Check for denormalized, exponent should be zero. Left *
350 * operand was normalized, so extent (guard, round) was zero */
351 if(!(underflowtrap = Is_underflowtrap_enabled()) &&
352 result_exponent==0) goto underflow;
353
354 /* Shift extension to complete one bit of normalization and
355 * update exponent. */
356 Ext_leftshiftby1(extent);
357
358 /* Discover first one bit to determine shift amount. Use a
359 * modified binary search. We have already shifted the result
360 * one position right and still not found a one so the remainder
361 * of the extension must be zero and simplifies rounding. */
362 /* Scan bytes */
363 while(Dbl_iszero_hiddenhigh7mantissa(resultp1))
364 {
365 Dbl_leftshiftby8(resultp1,resultp2);
366 if((result_exponent -= 8) <= 0 && !underflowtrap)
367 goto underflow;
368 }
369 /* Now narrow it down to the nibble */
370 if(Dbl_iszero_hiddenhigh3mantissa(resultp1))
371 {
372 /* The lower nibble contains the normalizing one */
373 Dbl_leftshiftby4(resultp1,resultp2);
374 if((result_exponent -= 4) <= 0 && !underflowtrap)
375 goto underflow;
376 }
377 /* Select case were first bit is set (already normalized)
378 * otherwise select the proper shift. */
379 if((jumpsize = Dbl_hiddenhigh3mantissa(resultp1)) > 7)
380 {
381 /* Already normalized */
382 if(result_exponent <= 0) goto underflow;
383 Dbl_set_sign(resultp1,/*using*/sign_save);
384 Dbl_set_exponent(resultp1,/*using*/result_exponent);
385 Dbl_copytoptr(resultp1,resultp2,dstptr);
386 return(NOEXCEPTION);
387 }
388 Dbl_sethigh4bits(resultp1,/*using*/sign_save);
389 switch(jumpsize)
390 {
391 case 1:
392 {
393 Dbl_leftshiftby3(resultp1,resultp2);
394 result_exponent -= 3;
395 break;
396 }
397 case 2:
398 case 3:
399 {
400 Dbl_leftshiftby2(resultp1,resultp2);
401 result_exponent -= 2;
402 break;
403 }
404 case 4:
405 case 5:
406 case 6:
407 case 7:
408 {
409 Dbl_leftshiftby1(resultp1,resultp2);
410 result_exponent -= 1;
411 break;
412 }
413 }
414 if(result_exponent > 0)
415 {
416 Dbl_set_exponent(resultp1,/*using*/result_exponent);
417 Dbl_copytoptr(resultp1,resultp2,dstptr);
418 return(NOEXCEPTION); /* Sign bit is already set */
419 }
420 /* Fixup potential underflows */
421 underflow:
422 if(Is_underflowtrap_enabled())
423 {
424 Dbl_set_sign(resultp1,sign_save);
425 Dbl_setwrapped_exponent(resultp1,result_exponent,unfl);
426 Dbl_copytoptr(resultp1,resultp2,dstptr);
427 /* inexact = false */
428 return(UNDERFLOWEXCEPTION);
429 }
430 /*
431 * Since we cannot get an inexact denormalized result,
432 * we can now return.
433 */
434 Dbl_fix_overshift(resultp1,resultp2,(1-result_exponent),extent);
435 Dbl_clear_signexponent(resultp1);
436 Dbl_set_sign(resultp1,sign_save);
437 Dbl_copytoptr(resultp1,resultp2,dstptr);
438 return(NOEXCEPTION);
439 } /* end if(hidden...)... */
440 /* Fall through and round */
441 } /* end if(save >= 0)... */
442 else
443 {
444 /* Subtract magnitudes */
445 Dbl_addition(leftp1,leftp2,rightp1,rightp2,/*to*/resultp1,resultp2);
446 if(Dbl_isone_hiddenoverflow(resultp1))
447 {
448 /* Prenormalization required. */
449 Dbl_rightshiftby1_withextent(resultp2,extent,extent);
450 Dbl_arithrightshiftby1(resultp1,resultp2);
451 result_exponent++;
452 } /* end if hiddenoverflow... */
453 } /* end else ...subtract magnitudes... */
454
455 /* Round the result. If the extension is all zeros,then the result is
456 * exact. Otherwise round in the correct direction. No underflow is
457 * possible. If a postnormalization is necessary, then the mantissa is
458 * all zeros so no shift is needed. */
459 round:
460 if(Ext_isnotzero(extent))
461 {
462 inexact = true;
463 switch(Rounding_mode())
464 {
465 case ROUNDNEAREST: /* The default. */
466 if(Ext_isone_sign(extent))
467 {
468 /* at least 1/2 ulp */
469 if(Ext_isnotzero_lower(extent) ||
470 Dbl_isone_lowmantissap2(resultp2))
471 {
472 /* either exactly half way and odd or more than 1/2ulp */
473 Dbl_increment(resultp1,resultp2);
474 }
475 }
476 break;
477
478 case ROUNDPLUS:
479 if(Dbl_iszero_sign(resultp1))
480 {
481 /* Round up positive results */
482 Dbl_increment(resultp1,resultp2);
483 }
484 break;
485
486 case ROUNDMINUS:
487 if(Dbl_isone_sign(resultp1))
488 {
489 /* Round down negative results */
490 Dbl_increment(resultp1,resultp2);
491 }
492
493 case ROUNDZERO:;
494 /* truncate is simple */
495 } /* end switch... */
496 if(Dbl_isone_hiddenoverflow(resultp1)) result_exponent++;
497 }
498 if(result_exponent == DBL_INFINITY_EXPONENT)
499 {
500 /* Overflow */
501 if(Is_overflowtrap_enabled())
502 {
503 Dbl_setwrapped_exponent(resultp1,result_exponent,ovfl);
504 Dbl_copytoptr(resultp1,resultp2,dstptr);
505 if (inexact) {
506 if (Is_inexacttrap_enabled())
507 return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
508 else
509 Set_inexactflag();
510 }
511 return(OVERFLOWEXCEPTION);
512 }
513 else
514 {
515 inexact = true;
516 Set_overflowflag();
517 Dbl_setoverflow(resultp1,resultp2);
518 }
519 }
520 else Dbl_set_exponent(resultp1,result_exponent);
521 Dbl_copytoptr(resultp1,resultp2,dstptr);
522 if(inexact) {
523 if(Is_inexacttrap_enabled())
524 return(INEXACTEXCEPTION);
525 else
526 Set_inexactflag();
527 }
528 return(NOEXCEPTION);
529 }
530