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