?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Linux 5.8-rc4
[linux/fpc-iii.git] / arch / parisc / math-emu / fmpyfadd.c
blob029f9bbfe7b5b93145280f8ad6e06669d1c420e3
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/fmpyfadd.c $Revision: 1.1 $
13 *
14 * Purpose:
15 * Double Floating-point Multiply Fused Add
16 * Double Floating-point Multiply Negate Fused Add
17 * Single Floating-point Multiply Fused Add
18 * Single Floating-point Multiply Negate Fused Add
19 *
20 * External Interfaces:
21 * dbl_fmpyfadd(src1ptr,src2ptr,src3ptr,status,dstptr)
22 * dbl_fmpynfadd(src1ptr,src2ptr,src3ptr,status,dstptr)
23 * sgl_fmpyfadd(src1ptr,src2ptr,src3ptr,status,dstptr)
24 * sgl_fmpynfadd(src1ptr,src2ptr,src3ptr,status,dstptr)
25 *
26 * Internal Interfaces:
27 *
28 * Theory:
29 * <<please update with a overview of the operation of this file>>
30 *
31 * END_DESC
32 */
33
34
35 #include "float.h"
36 #include "sgl_float.h"
37 #include "dbl_float.h"
38
39
40 /*
41 * Double Floating-point Multiply Fused Add
42 */
43
44 int
45 dbl_fmpyfadd(
46 dbl_floating_point *src1ptr,
47 dbl_floating_point *src2ptr,
48 dbl_floating_point *src3ptr,
49 unsigned int *status,
50 dbl_floating_point *dstptr)
51 {
52 unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2, opnd3p1, opnd3p2;
53 register unsigned int tmpresp1, tmpresp2, tmpresp3, tmpresp4;
54 unsigned int rightp1, rightp2, rightp3, rightp4;
55 unsigned int resultp1, resultp2 = 0, resultp3 = 0, resultp4 = 0;
56 register int mpy_exponent, add_exponent, count;
57 boolean inexact = FALSE, is_tiny = FALSE;
58
59 unsigned int signlessleft1, signlessright1, save;
60 register int result_exponent, diff_exponent;
61 int sign_save, jumpsize;
62
63 Dbl_copyfromptr(src1ptr,opnd1p1,opnd1p2);
64 Dbl_copyfromptr(src2ptr,opnd2p1,opnd2p2);
65 Dbl_copyfromptr(src3ptr,opnd3p1,opnd3p2);
66
67 /*
68 * set sign bit of result of multiply
69 */
70 if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1))
71 Dbl_setnegativezerop1(resultp1);
72 else Dbl_setzerop1(resultp1);
73
74 /*
75 * Generate multiply exponent
76 */
77 mpy_exponent = Dbl_exponent(opnd1p1) + Dbl_exponent(opnd2p1) - DBL_BIAS;
78
79 /*
80 * check first operand for NaN's or infinity
81 */
82 if (Dbl_isinfinity_exponent(opnd1p1)) {
83 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
84 if (Dbl_isnotnan(opnd2p1,opnd2p2) &&
85 Dbl_isnotnan(opnd3p1,opnd3p2)) {
86 if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
87 /*
88 * invalid since operands are infinity
89 * and zero
90 */
91 if (Is_invalidtrap_enabled())
92 return(OPC_2E_INVALIDEXCEPTION);
93 Set_invalidflag();
94 Dbl_makequietnan(resultp1,resultp2);
95 Dbl_copytoptr(resultp1,resultp2,dstptr);
96 return(NOEXCEPTION);
97 }
98 /*
99 * Check third operand for infinity with a
100 * sign opposite of the multiply result
101 */
102 if (Dbl_isinfinity(opnd3p1,opnd3p2) &&
103 (Dbl_sign(resultp1) ^ Dbl_sign(opnd3p1))) {
104 /*
105 * invalid since attempting a magnitude
106 * subtraction of infinities
107 */
108 if (Is_invalidtrap_enabled())
109 return(OPC_2E_INVALIDEXCEPTION);
110 Set_invalidflag();
111 Dbl_makequietnan(resultp1,resultp2);
112 Dbl_copytoptr(resultp1,resultp2,dstptr);
113 return(NOEXCEPTION);
114 }
115
116 /*
117 * return infinity
118 */
119 Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
120 Dbl_copytoptr(resultp1,resultp2,dstptr);
121 return(NOEXCEPTION);
122 }
123 }
124 else {
125 /*
126 * is NaN; signaling or quiet?
127 */
128 if (Dbl_isone_signaling(opnd1p1)) {
129 /* trap if INVALIDTRAP enabled */
130 if (Is_invalidtrap_enabled())
131 return(OPC_2E_INVALIDEXCEPTION);
132 /* make NaN quiet */
133 Set_invalidflag();
134 Dbl_set_quiet(opnd1p1);
135 }
136 /*
137 * is second operand a signaling NaN?
138 */
139 else if (Dbl_is_signalingnan(opnd2p1)) {
140 /* trap if INVALIDTRAP enabled */
141 if (Is_invalidtrap_enabled())
142 return(OPC_2E_INVALIDEXCEPTION);
143 /* make NaN quiet */
144 Set_invalidflag();
145 Dbl_set_quiet(opnd2p1);
146 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
147 return(NOEXCEPTION);
148 }
149 /*
150 * is third operand a signaling NaN?
151 */
152 else if (Dbl_is_signalingnan(opnd3p1)) {
153 /* trap if INVALIDTRAP enabled */
154 if (Is_invalidtrap_enabled())
155 return(OPC_2E_INVALIDEXCEPTION);
156 /* make NaN quiet */
157 Set_invalidflag();
158 Dbl_set_quiet(opnd3p1);
159 Dbl_copytoptr(opnd3p1,opnd3p2,dstptr);
160 return(NOEXCEPTION);
161 }
162 /*
163 * return quiet NaN
164 */
165 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
166 return(NOEXCEPTION);
167 }
168 }
169
170 /*
171 * check second operand for NaN's or infinity
172 */
173 if (Dbl_isinfinity_exponent(opnd2p1)) {
174 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
175 if (Dbl_isnotnan(opnd3p1,opnd3p2)) {
176 if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) {
177 /*
178 * invalid since multiply operands are
179 * zero & infinity
180 */
181 if (Is_invalidtrap_enabled())
182 return(OPC_2E_INVALIDEXCEPTION);
183 Set_invalidflag();
184 Dbl_makequietnan(opnd2p1,opnd2p2);
185 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
186 return(NOEXCEPTION);
187 }
188
189 /*
190 * Check third operand for infinity with a
191 * sign opposite of the multiply result
192 */
193 if (Dbl_isinfinity(opnd3p1,opnd3p2) &&
194 (Dbl_sign(resultp1) ^ Dbl_sign(opnd3p1))) {
195 /*
196 * invalid since attempting a magnitude
197 * subtraction of infinities
198 */
199 if (Is_invalidtrap_enabled())
200 return(OPC_2E_INVALIDEXCEPTION);
201 Set_invalidflag();
202 Dbl_makequietnan(resultp1,resultp2);
203 Dbl_copytoptr(resultp1,resultp2,dstptr);
204 return(NOEXCEPTION);
205 }
206
207 /*
208 * return infinity
209 */
210 Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
211 Dbl_copytoptr(resultp1,resultp2,dstptr);
212 return(NOEXCEPTION);
213 }
214 }
215 else {
216 /*
217 * is NaN; signaling or quiet?
218 */
219 if (Dbl_isone_signaling(opnd2p1)) {
220 /* trap if INVALIDTRAP enabled */
221 if (Is_invalidtrap_enabled())
222 return(OPC_2E_INVALIDEXCEPTION);
223 /* make NaN quiet */
224 Set_invalidflag();
225 Dbl_set_quiet(opnd2p1);
226 }
227 /*
228 * is third operand a signaling NaN?
229 */
230 else if (Dbl_is_signalingnan(opnd3p1)) {
231 /* trap if INVALIDTRAP enabled */
232 if (Is_invalidtrap_enabled())
233 return(OPC_2E_INVALIDEXCEPTION);
234 /* make NaN quiet */
235 Set_invalidflag();
236 Dbl_set_quiet(opnd3p1);
237 Dbl_copytoptr(opnd3p1,opnd3p2,dstptr);
238 return(NOEXCEPTION);
239 }
240 /*
241 * return quiet NaN
242 */
243 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
244 return(NOEXCEPTION);
245 }
246 }
247
248 /*
249 * check third operand for NaN's or infinity
250 */
251 if (Dbl_isinfinity_exponent(opnd3p1)) {
252 if (Dbl_iszero_mantissa(opnd3p1,opnd3p2)) {
253 /* return infinity */
254 Dbl_copytoptr(opnd3p1,opnd3p2,dstptr);
255 return(NOEXCEPTION);
256 } else {
257 /*
258 * is NaN; signaling or quiet?
259 */
260 if (Dbl_isone_signaling(opnd3p1)) {
261 /* trap if INVALIDTRAP enabled */
262 if (Is_invalidtrap_enabled())
263 return(OPC_2E_INVALIDEXCEPTION);
264 /* make NaN quiet */
265 Set_invalidflag();
266 Dbl_set_quiet(opnd3p1);
267 }
268 /*
269 * return quiet NaN
270 */
271 Dbl_copytoptr(opnd3p1,opnd3p2,dstptr);
272 return(NOEXCEPTION);
273 }
274 }
275
276 /*
277 * Generate multiply mantissa
278 */
279 if (Dbl_isnotzero_exponent(opnd1p1)) {
280 /* set hidden bit */
281 Dbl_clear_signexponent_set_hidden(opnd1p1);
282 }
283 else {
284 /* check for zero */
285 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
286 /*
287 * Perform the add opnd3 with zero here.
288 */
289 if (Dbl_iszero_exponentmantissa(opnd3p1,opnd3p2)) {
290 if (Is_rounding_mode(ROUNDMINUS)) {
291 Dbl_or_signs(opnd3p1,resultp1);
292 } else {
293 Dbl_and_signs(opnd3p1,resultp1);
294 }
295 }
296 /*
297 * Now let's check for trapped underflow case.
298 */
299 else if (Dbl_iszero_exponent(opnd3p1) &&
300 Is_underflowtrap_enabled()) {
301 /* need to normalize results mantissa */
302 sign_save = Dbl_signextendedsign(opnd3p1);
303 result_exponent = 0;
304 Dbl_leftshiftby1(opnd3p1,opnd3p2);
305 Dbl_normalize(opnd3p1,opnd3p2,result_exponent);
306 Dbl_set_sign(opnd3p1,/*using*/sign_save);
307 Dbl_setwrapped_exponent(opnd3p1,result_exponent,
308 unfl);
309 Dbl_copytoptr(opnd3p1,opnd3p2,dstptr);
310 /* inexact = FALSE */
311 return(OPC_2E_UNDERFLOWEXCEPTION);
312 }
313 Dbl_copytoptr(opnd3p1,opnd3p2,dstptr);
314 return(NOEXCEPTION);
315 }
316 /* is denormalized, adjust exponent */
317 Dbl_clear_signexponent(opnd1p1);
318 Dbl_leftshiftby1(opnd1p1,opnd1p2);
319 Dbl_normalize(opnd1p1,opnd1p2,mpy_exponent);
320 }
321 /* opnd2 needs to have hidden bit set with msb in hidden bit */
322 if (Dbl_isnotzero_exponent(opnd2p1)) {
323 Dbl_clear_signexponent_set_hidden(opnd2p1);
324 }
325 else {
326 /* check for zero */
327 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
328 /*
329 * Perform the add opnd3 with zero here.
330 */
331 if (Dbl_iszero_exponentmantissa(opnd3p1,opnd3p2)) {
332 if (Is_rounding_mode(ROUNDMINUS)) {
333 Dbl_or_signs(opnd3p1,resultp1);
334 } else {
335 Dbl_and_signs(opnd3p1,resultp1);
336 }
337 }
338 /*
339 * Now let's check for trapped underflow case.
340 */
341 else if (Dbl_iszero_exponent(opnd3p1) &&
342 Is_underflowtrap_enabled()) {
343 /* need to normalize results mantissa */
344 sign_save = Dbl_signextendedsign(opnd3p1);
345 result_exponent = 0;
346 Dbl_leftshiftby1(opnd3p1,opnd3p2);
347 Dbl_normalize(opnd3p1,opnd3p2,result_exponent);
348 Dbl_set_sign(opnd3p1,/*using*/sign_save);
349 Dbl_setwrapped_exponent(opnd3p1,result_exponent,
350 unfl);
351 Dbl_copytoptr(opnd3p1,opnd3p2,dstptr);
352 /* inexact = FALSE */
353 return(OPC_2E_UNDERFLOWEXCEPTION);
354 }
355 Dbl_copytoptr(opnd3p1,opnd3p2,dstptr);
356 return(NOEXCEPTION);
357 }
358 /* is denormalized; want to normalize */
359 Dbl_clear_signexponent(opnd2p1);
360 Dbl_leftshiftby1(opnd2p1,opnd2p2);
361 Dbl_normalize(opnd2p1,opnd2p2,mpy_exponent);
362 }
363
364 /* Multiply the first two source mantissas together */
365
366 /*
367 * The intermediate result will be kept in tmpres,
368 * which needs enough room for 106 bits of mantissa,
369 * so lets call it a Double extended.
370 */
371 Dblext_setzero(tmpresp1,tmpresp2,tmpresp3,tmpresp4);
372
373 /*
374 * Four bits at a time are inspected in each loop, and a
375 * simple shift and add multiply algorithm is used.
376 */
377 for (count = DBL_P-1; count >= 0; count -= 4) {
378 Dblext_rightshiftby4(tmpresp1,tmpresp2,tmpresp3,tmpresp4);
379 if (Dbit28p2(opnd1p2)) {
380 /* Fourword_add should be an ADD followed by 3 ADDC's */
381 Fourword_add(tmpresp1, tmpresp2, tmpresp3, tmpresp4,
382 opnd2p1<<3 | opnd2p2>>29, opnd2p2<<3, 0, 0);
383 }
384 if (Dbit29p2(opnd1p2)) {
385 Fourword_add(tmpresp1, tmpresp2, tmpresp3, tmpresp4,
386 opnd2p1<<2 | opnd2p2>>30, opnd2p2<<2, 0, 0);
387 }
388 if (Dbit30p2(opnd1p2)) {
389 Fourword_add(tmpresp1, tmpresp2, tmpresp3, tmpresp4,
390 opnd2p1<<1 | opnd2p2>>31, opnd2p2<<1, 0, 0);
391 }
392 if (Dbit31p2(opnd1p2)) {
393 Fourword_add(tmpresp1, tmpresp2, tmpresp3, tmpresp4,
394 opnd2p1, opnd2p2, 0, 0);
395 }
396 Dbl_rightshiftby4(opnd1p1,opnd1p2);
397 }
398 if (Is_dexthiddenoverflow(tmpresp1)) {
399 /* result mantissa >= 2 (mantissa overflow) */
400 mpy_exponent++;
401 Dblext_rightshiftby1(tmpresp1,tmpresp2,tmpresp3,tmpresp4);
402 }
403
404 /*
405 * Restore the sign of the mpy result which was saved in resultp1.
406 * The exponent will continue to be kept in mpy_exponent.
407 */
408 Dblext_set_sign(tmpresp1,Dbl_sign(resultp1));
409
410 /*
411 * No rounding is required, since the result of the multiply
412 * is exact in the extended format.
413 */
414
415 /*
416 * Now we are ready to perform the add portion of the operation.
417 *
418 * The exponents need to be kept as integers for now, since the
419 * multiply result might not fit into the exponent field. We
420 * can't overflow or underflow because of this yet, since the
421 * add could bring the final result back into range.
422 */
423 add_exponent = Dbl_exponent(opnd3p1);
424
425 /*
426 * Check for denormalized or zero add operand.
427 */
428 if (add_exponent == 0) {
429 /* check for zero */
430 if (Dbl_iszero_mantissa(opnd3p1,opnd3p2)) {
431 /* right is zero */
432 /* Left can't be zero and must be result.
433 *
434 * The final result is now in tmpres and mpy_exponent,
435 * and needs to be rounded and squeezed back into
436 * double precision format from double extended.
437 */
438 result_exponent = mpy_exponent;
439 Dblext_copy(tmpresp1,tmpresp2,tmpresp3,tmpresp4,
440 resultp1,resultp2,resultp3,resultp4);
441 sign_save = Dbl_signextendedsign(resultp1);/*save sign*/
442 goto round;
443 }
444
445 /*
446 * Neither are zeroes.
447 * Adjust exponent and normalize add operand.
448 */
449 sign_save = Dbl_signextendedsign(opnd3p1); /* save sign */
450 Dbl_clear_signexponent(opnd3p1);
451 Dbl_leftshiftby1(opnd3p1,opnd3p2);
452 Dbl_normalize(opnd3p1,opnd3p2,add_exponent);
453 Dbl_set_sign(opnd3p1,sign_save); /* restore sign */
454 } else {
455 Dbl_clear_exponent_set_hidden(opnd3p1);
456 }
457 /*
458 * Copy opnd3 to the double extended variable called right.
459 */
460 Dbl_copyto_dblext(opnd3p1,opnd3p2,rightp1,rightp2,rightp3,rightp4);
461
462 /*
463 * A zero "save" helps discover equal operands (for later),
464 * and is used in swapping operands (if needed).
465 */
466 Dblext_xortointp1(tmpresp1,rightp1,/*to*/save);
467
468 /*
469 * Compare magnitude of operands.
470 */
471 Dblext_copytoint_exponentmantissap1(tmpresp1,signlessleft1);
472 Dblext_copytoint_exponentmantissap1(rightp1,signlessright1);
473 if (mpy_exponent < add_exponent || mpy_exponent == add_exponent &&
474 Dblext_ismagnitudeless(tmpresp2,rightp2,signlessleft1,signlessright1)){
475 /*
476 * Set the left operand to the larger one by XOR swap.
477 * First finish the first word "save".
478 */
479 Dblext_xorfromintp1(save,rightp1,/*to*/rightp1);
480 Dblext_xorfromintp1(save,tmpresp1,/*to*/tmpresp1);
481 Dblext_swap_lower(tmpresp2,tmpresp3,tmpresp4,
482 rightp2,rightp3,rightp4);
483 /* also setup exponents used in rest of routine */
484 diff_exponent = add_exponent - mpy_exponent;
485 result_exponent = add_exponent;
486 } else {
487 /* also setup exponents used in rest of routine */
488 diff_exponent = mpy_exponent - add_exponent;
489 result_exponent = mpy_exponent;
490 }
491 /* Invariant: left is not smaller than right. */
492
493 /*
494 * Special case alignment of operands that would force alignment
495 * beyond the extent of the extension. A further optimization
496 * could special case this but only reduces the path length for
497 * this infrequent case.
498 */
499 if (diff_exponent > DBLEXT_THRESHOLD) {
500 diff_exponent = DBLEXT_THRESHOLD;
501 }
502
503 /* Align right operand by shifting it to the right */
504 Dblext_clear_sign(rightp1);
505 Dblext_right_align(rightp1,rightp2,rightp3,rightp4,
506 /*shifted by*/diff_exponent);
507
508 /* Treat sum and difference of the operands separately. */
509 if ((int)save < 0) {
510 /*
511 * Difference of the two operands. Overflow can occur if the
512 * multiply overflowed. A borrow can occur out of the hidden
513 * bit and force a post normalization phase.
514 */
515 Dblext_subtract(tmpresp1,tmpresp2,tmpresp3,tmpresp4,
516 rightp1,rightp2,rightp3,rightp4,
517 resultp1,resultp2,resultp3,resultp4);
518 sign_save = Dbl_signextendedsign(resultp1);
519 if (Dbl_iszero_hidden(resultp1)) {
520 /* Handle normalization */
521 /* A straightforward algorithm would now shift the
522 * result and extension left until the hidden bit
523 * becomes one. Not all of the extension bits need
524 * participate in the shift. Only the two most
525 * significant bits (round and guard) are needed.
526 * If only a single shift is needed then the guard
527 * bit becomes a significant low order bit and the
528 * extension must participate in the rounding.
529 * If more than a single shift is needed, then all
530 * bits to the right of the guard bit are zeros,
531 * and the guard bit may or may not be zero. */
532 Dblext_leftshiftby1(resultp1,resultp2,resultp3,
533 resultp4);
534
535 /* Need to check for a zero result. The sign and
536 * exponent fields have already been zeroed. The more
537 * efficient test of the full object can be used.
538 */
539 if(Dblext_iszero(resultp1,resultp2,resultp3,resultp4)){
540 /* Must have been "x-x" or "x+(-x)". */
541 if (Is_rounding_mode(ROUNDMINUS))
542 Dbl_setone_sign(resultp1);
543 Dbl_copytoptr(resultp1,resultp2,dstptr);
544 return(NOEXCEPTION);
545 }
546 result_exponent--;
547
548 /* Look to see if normalization is finished. */
549 if (Dbl_isone_hidden(resultp1)) {
550 /* No further normalization is needed */
551 goto round;
552 }
553
554 /* Discover first one bit to determine shift amount.
555 * Use a modified binary search. We have already
556 * shifted the result one position right and still
557 * not found a one so the remainder of the extension
558 * must be zero and simplifies rounding. */
559 /* Scan bytes */
560 while (Dbl_iszero_hiddenhigh7mantissa(resultp1)) {
561 Dblext_leftshiftby8(resultp1,resultp2,resultp3,resultp4);
562 result_exponent -= 8;
563 }
564 /* Now narrow it down to the nibble */
565 if (Dbl_iszero_hiddenhigh3mantissa(resultp1)) {
566 /* The lower nibble contains the
567 * normalizing one */
568 Dblext_leftshiftby4(resultp1,resultp2,resultp3,resultp4);
569 result_exponent -= 4;
570 }
571 /* Select case where first bit is set (already
572 * normalized) otherwise select the proper shift. */
573 jumpsize = Dbl_hiddenhigh3mantissa(resultp1);
574 if (jumpsize <= 7) switch(jumpsize) {
575 case 1:
576 Dblext_leftshiftby3(resultp1,resultp2,resultp3,
577 resultp4);
578 result_exponent -= 3;
579 break;
580 case 2:
581 case 3:
582 Dblext_leftshiftby2(resultp1,resultp2,resultp3,
583 resultp4);
584 result_exponent -= 2;
585 break;
586 case 4:
587 case 5:
588 case 6:
589 case 7:
590 Dblext_leftshiftby1(resultp1,resultp2,resultp3,
591 resultp4);
592 result_exponent -= 1;
593 break;
594 }
595 } /* end if (hidden...)... */
596 /* Fall through and round */
597 } /* end if (save < 0)... */
598 else {
599 /* Add magnitudes */
600 Dblext_addition(tmpresp1,tmpresp2,tmpresp3,tmpresp4,
601 rightp1,rightp2,rightp3,rightp4,
602 /*to*/resultp1,resultp2,resultp3,resultp4);
603 sign_save = Dbl_signextendedsign(resultp1);
604 if (Dbl_isone_hiddenoverflow(resultp1)) {
605 /* Prenormalization required. */
606 Dblext_arithrightshiftby1(resultp1,resultp2,resultp3,
607 resultp4);
608 result_exponent++;
609 } /* end if hiddenoverflow... */
610 } /* end else ...add magnitudes... */
611
612 /* Round the result. If the extension and lower two words are
613 * all zeros, then the result is exact. Otherwise round in the
614 * correct direction. Underflow is possible. If a postnormalization
615 * is necessary, then the mantissa is all zeros so no shift is needed.
616 */
617 round:
618 if (result_exponent <= 0 && !Is_underflowtrap_enabled()) {
619 Dblext_denormalize(resultp1,resultp2,resultp3,resultp4,
620 result_exponent,is_tiny);
621 }
622 Dbl_set_sign(resultp1,/*using*/sign_save);
623 if (Dblext_isnotzero_mantissap3(resultp3) ||
624 Dblext_isnotzero_mantissap4(resultp4)) {
625 inexact = TRUE;
626 switch(Rounding_mode()) {
627 case ROUNDNEAREST: /* The default. */
628 if (Dblext_isone_highp3(resultp3)) {
629 /* at least 1/2 ulp */
630 if (Dblext_isnotzero_low31p3(resultp3) ||
631 Dblext_isnotzero_mantissap4(resultp4) ||
632 Dblext_isone_lowp2(resultp2)) {
633 /* either exactly half way and odd or
634 * more than 1/2ulp */
635 Dbl_increment(resultp1,resultp2);
636 }
637 }
638 break;
639
640 case ROUNDPLUS:
641 if (Dbl_iszero_sign(resultp1)) {
642 /* Round up positive results */
643 Dbl_increment(resultp1,resultp2);
644 }
645 break;
646
647 case ROUNDMINUS:
648 if (Dbl_isone_sign(resultp1)) {
649 /* Round down negative results */
650 Dbl_increment(resultp1,resultp2);
651 }
652
653 case ROUNDZERO:;
654 /* truncate is simple */
655 } /* end switch... */
656 if (Dbl_isone_hiddenoverflow(resultp1)) result_exponent++;
657 }
658 if (result_exponent >= DBL_INFINITY_EXPONENT) {
659 /* trap if OVERFLOWTRAP enabled */
660 if (Is_overflowtrap_enabled()) {
661 /*
662 * Adjust bias of result
663 */
664 Dbl_setwrapped_exponent(resultp1,result_exponent,ovfl);
665 Dbl_copytoptr(resultp1,resultp2,dstptr);
666 if (inexact)
667 if (Is_inexacttrap_enabled())
668 return (OPC_2E_OVERFLOWEXCEPTION |
669 OPC_2E_INEXACTEXCEPTION);
670 else Set_inexactflag();
671 return (OPC_2E_OVERFLOWEXCEPTION);
672 }
673 inexact = TRUE;
674 Set_overflowflag();
675 /* set result to infinity or largest number */
676 Dbl_setoverflow(resultp1,resultp2);
677
678 } else if (result_exponent <= 0) { /* underflow case */
679 if (Is_underflowtrap_enabled()) {
680 /*
681 * Adjust bias of result
682 */
683 Dbl_setwrapped_exponent(resultp1,result_exponent,unfl);
684 Dbl_copytoptr(resultp1,resultp2,dstptr);
685 if (inexact)
686 if (Is_inexacttrap_enabled())
687 return (OPC_2E_UNDERFLOWEXCEPTION |
688 OPC_2E_INEXACTEXCEPTION);
689 else Set_inexactflag();
690 return(OPC_2E_UNDERFLOWEXCEPTION);
691 }
692 else if (inexact && is_tiny) Set_underflowflag();
693 }
694 else Dbl_set_exponent(resultp1,result_exponent);
695 Dbl_copytoptr(resultp1,resultp2,dstptr);
696 if (inexact)
697 if (Is_inexacttrap_enabled()) return(OPC_2E_INEXACTEXCEPTION);
698 else Set_inexactflag();
699 return(NOEXCEPTION);
700 }
701
702 /*
703 * Double Floating-point Multiply Negate Fused Add
704 */
705
706 dbl_fmpynfadd(src1ptr,src2ptr,src3ptr,status,dstptr)
707
708 dbl_floating_point *src1ptr, *src2ptr, *src3ptr, *dstptr;
709 unsigned int *status;
710 {
711 unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2, opnd3p1, opnd3p2;
712 register unsigned int tmpresp1, tmpresp2, tmpresp3, tmpresp4;
713 unsigned int rightp1, rightp2, rightp3, rightp4;
714 unsigned int resultp1, resultp2 = 0, resultp3 = 0, resultp4 = 0;
715 register int mpy_exponent, add_exponent, count;
716 boolean inexact = FALSE, is_tiny = FALSE;
717
718 unsigned int signlessleft1, signlessright1, save;
719 register int result_exponent, diff_exponent;
720 int sign_save, jumpsize;
721
722 Dbl_copyfromptr(src1ptr,opnd1p1,opnd1p2);
723 Dbl_copyfromptr(src2ptr,opnd2p1,opnd2p2);
724 Dbl_copyfromptr(src3ptr,opnd3p1,opnd3p2);
725
726 /*
727 * set sign bit of result of multiply
728 */
729 if (Dbl_sign(opnd1p1) ^ Dbl_sign(opnd2p1))
730 Dbl_setzerop1(resultp1);
731 else
732 Dbl_setnegativezerop1(resultp1);
733
734 /*
735 * Generate multiply exponent
736 */
737 mpy_exponent = Dbl_exponent(opnd1p1) + Dbl_exponent(opnd2p1) - DBL_BIAS;
738
739 /*
740 * check first operand for NaN's or infinity
741 */
742 if (Dbl_isinfinity_exponent(opnd1p1)) {
743 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
744 if (Dbl_isnotnan(opnd2p1,opnd2p2) &&
745 Dbl_isnotnan(opnd3p1,opnd3p2)) {
746 if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
747 /*
748 * invalid since operands are infinity
749 * and zero
750 */
751 if (Is_invalidtrap_enabled())
752 return(OPC_2E_INVALIDEXCEPTION);
753 Set_invalidflag();
754 Dbl_makequietnan(resultp1,resultp2);
755 Dbl_copytoptr(resultp1,resultp2,dstptr);
756 return(NOEXCEPTION);
757 }
758 /*
759 * Check third operand for infinity with a
760 * sign opposite of the multiply result
761 */
762 if (Dbl_isinfinity(opnd3p1,opnd3p2) &&
763 (Dbl_sign(resultp1) ^ Dbl_sign(opnd3p1))) {
764 /*
765 * invalid since attempting a magnitude
766 * subtraction of infinities
767 */
768 if (Is_invalidtrap_enabled())
769 return(OPC_2E_INVALIDEXCEPTION);
770 Set_invalidflag();
771 Dbl_makequietnan(resultp1,resultp2);
772 Dbl_copytoptr(resultp1,resultp2,dstptr);
773 return(NOEXCEPTION);
774 }
775
776 /*
777 * return infinity
778 */
779 Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
780 Dbl_copytoptr(resultp1,resultp2,dstptr);
781 return(NOEXCEPTION);
782 }
783 }
784 else {
785 /*
786 * is NaN; signaling or quiet?
787 */
788 if (Dbl_isone_signaling(opnd1p1)) {
789 /* trap if INVALIDTRAP enabled */
790 if (Is_invalidtrap_enabled())
791 return(OPC_2E_INVALIDEXCEPTION);
792 /* make NaN quiet */
793 Set_invalidflag();
794 Dbl_set_quiet(opnd1p1);
795 }
796 /*
797 * is second operand a signaling NaN?
798 */
799 else if (Dbl_is_signalingnan(opnd2p1)) {
800 /* trap if INVALIDTRAP enabled */
801 if (Is_invalidtrap_enabled())
802 return(OPC_2E_INVALIDEXCEPTION);
803 /* make NaN quiet */
804 Set_invalidflag();
805 Dbl_set_quiet(opnd2p1);
806 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
807 return(NOEXCEPTION);
808 }
809 /*
810 * is third operand a signaling NaN?
811 */
812 else if (Dbl_is_signalingnan(opnd3p1)) {
813 /* trap if INVALIDTRAP enabled */
814 if (Is_invalidtrap_enabled())
815 return(OPC_2E_INVALIDEXCEPTION);
816 /* make NaN quiet */
817 Set_invalidflag();
818 Dbl_set_quiet(opnd3p1);
819 Dbl_copytoptr(opnd3p1,opnd3p2,dstptr);
820 return(NOEXCEPTION);
821 }
822 /*
823 * return quiet NaN
824 */
825 Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
826 return(NOEXCEPTION);
827 }
828 }
829
830 /*
831 * check second operand for NaN's or infinity
832 */
833 if (Dbl_isinfinity_exponent(opnd2p1)) {
834 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
835 if (Dbl_isnotnan(opnd3p1,opnd3p2)) {
836 if (Dbl_iszero_exponentmantissa(opnd1p1,opnd1p2)) {
837 /*
838 * invalid since multiply operands are
839 * zero & infinity
840 */
841 if (Is_invalidtrap_enabled())
842 return(OPC_2E_INVALIDEXCEPTION);
843 Set_invalidflag();
844 Dbl_makequietnan(opnd2p1,opnd2p2);
845 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
846 return(NOEXCEPTION);
847 }
848
849 /*
850 * Check third operand for infinity with a
851 * sign opposite of the multiply result
852 */
853 if (Dbl_isinfinity(opnd3p1,opnd3p2) &&
854 (Dbl_sign(resultp1) ^ Dbl_sign(opnd3p1))) {
855 /*
856 * invalid since attempting a magnitude
857 * subtraction of infinities
858 */
859 if (Is_invalidtrap_enabled())
860 return(OPC_2E_INVALIDEXCEPTION);
861 Set_invalidflag();
862 Dbl_makequietnan(resultp1,resultp2);
863 Dbl_copytoptr(resultp1,resultp2,dstptr);
864 return(NOEXCEPTION);
865 }
866
867 /*
868 * return infinity
869 */
870 Dbl_setinfinity_exponentmantissa(resultp1,resultp2);
871 Dbl_copytoptr(resultp1,resultp2,dstptr);
872 return(NOEXCEPTION);
873 }
874 }
875 else {
876 /*
877 * is NaN; signaling or quiet?
878 */
879 if (Dbl_isone_signaling(opnd2p1)) {
880 /* trap if INVALIDTRAP enabled */
881 if (Is_invalidtrap_enabled())
882 return(OPC_2E_INVALIDEXCEPTION);
883 /* make NaN quiet */
884 Set_invalidflag();
885 Dbl_set_quiet(opnd2p1);
886 }
887 /*
888 * is third operand a signaling NaN?
889 */
890 else if (Dbl_is_signalingnan(opnd3p1)) {
891 /* trap if INVALIDTRAP enabled */
892 if (Is_invalidtrap_enabled())
893 return(OPC_2E_INVALIDEXCEPTION);
894 /* make NaN quiet */
895 Set_invalidflag();
896 Dbl_set_quiet(opnd3p1);
897 Dbl_copytoptr(opnd3p1,opnd3p2,dstptr);
898 return(NOEXCEPTION);
899 }
900 /*
901 * return quiet NaN
902 */
903 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
904 return(NOEXCEPTION);
905 }
906 }
907
908 /*
909 * check third operand for NaN's or infinity
910 */
911 if (Dbl_isinfinity_exponent(opnd3p1)) {
912 if (Dbl_iszero_mantissa(opnd3p1,opnd3p2)) {
913 /* return infinity */
914 Dbl_copytoptr(opnd3p1,opnd3p2,dstptr);
915 return(NOEXCEPTION);
916 } else {
917 /*
918 * is NaN; signaling or quiet?
919 */
920 if (Dbl_isone_signaling(opnd3p1)) {
921 /* trap if INVALIDTRAP enabled */
922 if (Is_invalidtrap_enabled())
923 return(OPC_2E_INVALIDEXCEPTION);
924 /* make NaN quiet */
925 Set_invalidflag();
926 Dbl_set_quiet(opnd3p1);
927 }
928 /*
929 * return quiet NaN
930 */
931 Dbl_copytoptr(opnd3p1,opnd3p2,dstptr);
932 return(NOEXCEPTION);
933 }
934 }
935
936 /*
937 * Generate multiply mantissa
938 */
939 if (Dbl_isnotzero_exponent(opnd1p1)) {
940 /* set hidden bit */
941 Dbl_clear_signexponent_set_hidden(opnd1p1);
942 }
943 else {
944 /* check for zero */
945 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
946 /*
947 * Perform the add opnd3 with zero here.
948 */
949 if (Dbl_iszero_exponentmantissa(opnd3p1,opnd3p2)) {
950 if (Is_rounding_mode(ROUNDMINUS)) {
951 Dbl_or_signs(opnd3p1,resultp1);
952 } else {
953 Dbl_and_signs(opnd3p1,resultp1);
954 }
955 }
956 /*
957 * Now let's check for trapped underflow case.
958 */
959 else if (Dbl_iszero_exponent(opnd3p1) &&
960 Is_underflowtrap_enabled()) {
961 /* need to normalize results mantissa */
962 sign_save = Dbl_signextendedsign(opnd3p1);
963 result_exponent = 0;
964 Dbl_leftshiftby1(opnd3p1,opnd3p2);
965 Dbl_normalize(opnd3p1,opnd3p2,result_exponent);
966 Dbl_set_sign(opnd3p1,/*using*/sign_save);
967 Dbl_setwrapped_exponent(opnd3p1,result_exponent,
968 unfl);
969 Dbl_copytoptr(opnd3p1,opnd3p2,dstptr);
970 /* inexact = FALSE */
971 return(OPC_2E_UNDERFLOWEXCEPTION);
972 }
973 Dbl_copytoptr(opnd3p1,opnd3p2,dstptr);
974 return(NOEXCEPTION);
975 }
976 /* is denormalized, adjust exponent */
977 Dbl_clear_signexponent(opnd1p1);
978 Dbl_leftshiftby1(opnd1p1,opnd1p2);
979 Dbl_normalize(opnd1p1,opnd1p2,mpy_exponent);
980 }
981 /* opnd2 needs to have hidden bit set with msb in hidden bit */
982 if (Dbl_isnotzero_exponent(opnd2p1)) {
983 Dbl_clear_signexponent_set_hidden(opnd2p1);
984 }
985 else {
986 /* check for zero */
987 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
988 /*
989 * Perform the add opnd3 with zero here.
990 */
991 if (Dbl_iszero_exponentmantissa(opnd3p1,opnd3p2)) {
992 if (Is_rounding_mode(ROUNDMINUS)) {
993 Dbl_or_signs(opnd3p1,resultp1);
994 } else {
995 Dbl_and_signs(opnd3p1,resultp1);
996 }
997 }
998 /*
999 * Now let's check for trapped underflow case.
1000 */
1001 else if (Dbl_iszero_exponent(opnd3p1) &&
1002 Is_underflowtrap_enabled()) {
1003 /* need to normalize results mantissa */
1004 sign_save = Dbl_signextendedsign(opnd3p1);
1005 result_exponent = 0;
1006 Dbl_leftshiftby1(opnd3p1,opnd3p2);
1007 Dbl_normalize(opnd3p1,opnd3p2,result_exponent);
1008 Dbl_set_sign(opnd3p1,/*using*/sign_save);
1009 Dbl_setwrapped_exponent(opnd3p1,result_exponent,
1010 unfl);
1011 Dbl_copytoptr(opnd3p1,opnd3p2,dstptr);
1012 /* inexact = FALSE */
1013 return(OPC_2E_UNDERFLOWEXCEPTION);
1014 }
1015 Dbl_copytoptr(opnd3p1,opnd3p2,dstptr);
1016 return(NOEXCEPTION);
1017 }
1018 /* is denormalized; want to normalize */
1019 Dbl_clear_signexponent(opnd2p1);
1020 Dbl_leftshiftby1(opnd2p1,opnd2p2);
1021 Dbl_normalize(opnd2p1,opnd2p2,mpy_exponent);
1022 }
1023
1024 /* Multiply the first two source mantissas together */
1025
1026 /*
1027 * The intermediate result will be kept in tmpres,
1028 * which needs enough room for 106 bits of mantissa,
1029 * so lets call it a Double extended.
1030 */
1031 Dblext_setzero(tmpresp1,tmpresp2,tmpresp3,tmpresp4);
1032
1033 /*
1034 * Four bits at a time are inspected in each loop, and a
1035 * simple shift and add multiply algorithm is used.
1036 */
1037 for (count = DBL_P-1; count >= 0; count -= 4) {
1038 Dblext_rightshiftby4(tmpresp1,tmpresp2,tmpresp3,tmpresp4);
1039 if (Dbit28p2(opnd1p2)) {
1040 /* Fourword_add should be an ADD followed by 3 ADDC's */
1041 Fourword_add(tmpresp1, tmpresp2, tmpresp3, tmpresp4,
1042 opnd2p1<<3 | opnd2p2>>29, opnd2p2<<3, 0, 0);
1043 }
1044 if (Dbit29p2(opnd1p2)) {
1045 Fourword_add(tmpresp1, tmpresp2, tmpresp3, tmpresp4,
1046 opnd2p1<<2 | opnd2p2>>30, opnd2p2<<2, 0, 0);
1047 }
1048 if (Dbit30p2(opnd1p2)) {
1049 Fourword_add(tmpresp1, tmpresp2, tmpresp3, tmpresp4,
1050 opnd2p1<<1 | opnd2p2>>31, opnd2p2<<1, 0, 0);
1051 }
1052 if (Dbit31p2(opnd1p2)) {
1053 Fourword_add(tmpresp1, tmpresp2, tmpresp3, tmpresp4,
1054 opnd2p1, opnd2p2, 0, 0);
1055 }
1056 Dbl_rightshiftby4(opnd1p1,opnd1p2);
1057 }
1058 if (Is_dexthiddenoverflow(tmpresp1)) {
1059 /* result mantissa >= 2 (mantissa overflow) */
1060 mpy_exponent++;
1061 Dblext_rightshiftby1(tmpresp1,tmpresp2,tmpresp3,tmpresp4);
1062 }
1063
1064 /*
1065 * Restore the sign of the mpy result which was saved in resultp1.
1066 * The exponent will continue to be kept in mpy_exponent.
1067 */
1068 Dblext_set_sign(tmpresp1,Dbl_sign(resultp1));
1069
1070 /*
1071 * No rounding is required, since the result of the multiply
1072 * is exact in the extended format.
1073 */
1074
1075 /*
1076 * Now we are ready to perform the add portion of the operation.
1077 *
1078 * The exponents need to be kept as integers for now, since the
1079 * multiply result might not fit into the exponent field. We
1080 * can't overflow or underflow because of this yet, since the
1081 * add could bring the final result back into range.
1082 */
1083 add_exponent = Dbl_exponent(opnd3p1);
1084
1085 /*
1086 * Check for denormalized or zero add operand.
1087 */
1088 if (add_exponent == 0) {
1089 /* check for zero */
1090 if (Dbl_iszero_mantissa(opnd3p1,opnd3p2)) {
1091 /* right is zero */
1092 /* Left can't be zero and must be result.
1093 *
1094 * The final result is now in tmpres and mpy_exponent,
1095 * and needs to be rounded and squeezed back into
1096 * double precision format from double extended.
1097 */
1098 result_exponent = mpy_exponent;
1099 Dblext_copy(tmpresp1,tmpresp2,tmpresp3,tmpresp4,
1100 resultp1,resultp2,resultp3,resultp4);
1101 sign_save = Dbl_signextendedsign(resultp1);/*save sign*/
1102 goto round;
1103 }
1104
1105 /*
1106 * Neither are zeroes.
1107 * Adjust exponent and normalize add operand.
1108 */
1109 sign_save = Dbl_signextendedsign(opnd3p1); /* save sign */
1110 Dbl_clear_signexponent(opnd3p1);
1111 Dbl_leftshiftby1(opnd3p1,opnd3p2);
1112 Dbl_normalize(opnd3p1,opnd3p2,add_exponent);
1113 Dbl_set_sign(opnd3p1,sign_save); /* restore sign */
1114 } else {
1115 Dbl_clear_exponent_set_hidden(opnd3p1);
1116 }
1117 /*
1118 * Copy opnd3 to the double extended variable called right.
1119 */
1120 Dbl_copyto_dblext(opnd3p1,opnd3p2,rightp1,rightp2,rightp3,rightp4);
1121
1122 /*
1123 * A zero "save" helps discover equal operands (for later),
1124 * and is used in swapping operands (if needed).
1125 */
1126 Dblext_xortointp1(tmpresp1,rightp1,/*to*/save);
1127
1128 /*
1129 * Compare magnitude of operands.
1130 */
1131 Dblext_copytoint_exponentmantissap1(tmpresp1,signlessleft1);
1132 Dblext_copytoint_exponentmantissap1(rightp1,signlessright1);
1133 if (mpy_exponent < add_exponent || mpy_exponent == add_exponent &&
1134 Dblext_ismagnitudeless(tmpresp2,rightp2,signlessleft1,signlessright1)){
1135 /*
1136 * Set the left operand to the larger one by XOR swap.
1137 * First finish the first word "save".
1138 */
1139 Dblext_xorfromintp1(save,rightp1,/*to*/rightp1);
1140 Dblext_xorfromintp1(save,tmpresp1,/*to*/tmpresp1);
1141 Dblext_swap_lower(tmpresp2,tmpresp3,tmpresp4,
1142 rightp2,rightp3,rightp4);
1143 /* also setup exponents used in rest of routine */
1144 diff_exponent = add_exponent - mpy_exponent;
1145 result_exponent = add_exponent;
1146 } else {
1147 /* also setup exponents used in rest of routine */
1148 diff_exponent = mpy_exponent - add_exponent;
1149 result_exponent = mpy_exponent;
1150 }
1151 /* Invariant: left is not smaller than right. */
1152
1153 /*
1154 * Special case alignment of operands that would force alignment
1155 * beyond the extent of the extension. A further optimization
1156 * could special case this but only reduces the path length for
1157 * this infrequent case.
1158 */
1159 if (diff_exponent > DBLEXT_THRESHOLD) {
1160 diff_exponent = DBLEXT_THRESHOLD;
1161 }
1162
1163 /* Align right operand by shifting it to the right */
1164 Dblext_clear_sign(rightp1);
1165 Dblext_right_align(rightp1,rightp2,rightp3,rightp4,
1166 /*shifted by*/diff_exponent);
1167
1168 /* Treat sum and difference of the operands separately. */
1169 if ((int)save < 0) {
1170 /*
1171 * Difference of the two operands. Overflow can occur if the
1172 * multiply overflowed. A borrow can occur out of the hidden
1173 * bit and force a post normalization phase.
1174 */
1175 Dblext_subtract(tmpresp1,tmpresp2,tmpresp3,tmpresp4,
1176 rightp1,rightp2,rightp3,rightp4,
1177 resultp1,resultp2,resultp3,resultp4);
1178 sign_save = Dbl_signextendedsign(resultp1);
1179 if (Dbl_iszero_hidden(resultp1)) {
1180 /* Handle normalization */
1181 /* A straightforward algorithm would now shift the
1182 * result and extension left until the hidden bit
1183 * becomes one. Not all of the extension bits need
1184 * participate in the shift. Only the two most
1185 * significant bits (round and guard) are needed.
1186 * If only a single shift is needed then the guard
1187 * bit becomes a significant low order bit and the
1188 * extension must participate in the rounding.
1189 * If more than a single shift is needed, then all
1190 * bits to the right of the guard bit are zeros,
1191 * and the guard bit may or may not be zero. */
1192 Dblext_leftshiftby1(resultp1,resultp2,resultp3,
1193 resultp4);
1194
1195 /* Need to check for a zero result. The sign and
1196 * exponent fields have already been zeroed. The more
1197 * efficient test of the full object can be used.
1198 */
1199 if (Dblext_iszero(resultp1,resultp2,resultp3,resultp4)) {
1200 /* Must have been "x-x" or "x+(-x)". */
1201 if (Is_rounding_mode(ROUNDMINUS))
1202 Dbl_setone_sign(resultp1);
1203 Dbl_copytoptr(resultp1,resultp2,dstptr);
1204 return(NOEXCEPTION);
1205 }
1206 result_exponent--;
1207
1208 /* Look to see if normalization is finished. */
1209 if (Dbl_isone_hidden(resultp1)) {
1210 /* No further normalization is needed */
1211 goto round;
1212 }
1213
1214 /* Discover first one bit to determine shift amount.
1215 * Use a modified binary search. We have already
1216 * shifted the result one position right and still
1217 * not found a one so the remainder of the extension
1218 * must be zero and simplifies rounding. */
1219 /* Scan bytes */
1220 while (Dbl_iszero_hiddenhigh7mantissa(resultp1)) {
1221 Dblext_leftshiftby8(resultp1,resultp2,resultp3,resultp4);
1222 result_exponent -= 8;
1223 }
1224 /* Now narrow it down to the nibble */
1225 if (Dbl_iszero_hiddenhigh3mantissa(resultp1)) {
1226 /* The lower nibble contains the
1227 * normalizing one */
1228 Dblext_leftshiftby4(resultp1,resultp2,resultp3,resultp4);
1229 result_exponent -= 4;
1230 }
1231 /* Select case where first bit is set (already
1232 * normalized) otherwise select the proper shift. */
1233 jumpsize = Dbl_hiddenhigh3mantissa(resultp1);
1234 if (jumpsize <= 7) switch(jumpsize) {
1235 case 1:
1236 Dblext_leftshiftby3(resultp1,resultp2,resultp3,
1237 resultp4);
1238 result_exponent -= 3;
1239 break;
1240 case 2:
1241 case 3:
1242 Dblext_leftshiftby2(resultp1,resultp2,resultp3,
1243 resultp4);
1244 result_exponent -= 2;
1245 break;
1246 case 4:
1247 case 5:
1248 case 6:
1249 case 7:
1250 Dblext_leftshiftby1(resultp1,resultp2,resultp3,
1251 resultp4);
1252 result_exponent -= 1;
1253 break;
1254 }
1255 } /* end if (hidden...)... */
1256 /* Fall through and round */
1257 } /* end if (save < 0)... */
1258 else {
1259 /* Add magnitudes */
1260 Dblext_addition(tmpresp1,tmpresp2,tmpresp3,tmpresp4,
1261 rightp1,rightp2,rightp3,rightp4,
1262 /*to*/resultp1,resultp2,resultp3,resultp4);
1263 sign_save = Dbl_signextendedsign(resultp1);
1264 if (Dbl_isone_hiddenoverflow(resultp1)) {
1265 /* Prenormalization required. */
1266 Dblext_arithrightshiftby1(resultp1,resultp2,resultp3,
1267 resultp4);
1268 result_exponent++;
1269 } /* end if hiddenoverflow... */
1270 } /* end else ...add magnitudes... */
1271
1272 /* Round the result. If the extension and lower two words are
1273 * all zeros, then the result is exact. Otherwise round in the
1274 * correct direction. Underflow is possible. If a postnormalization
1275 * is necessary, then the mantissa is all zeros so no shift is needed.
1276 */
1277 round:
1278 if (result_exponent <= 0 && !Is_underflowtrap_enabled()) {
1279 Dblext_denormalize(resultp1,resultp2,resultp3,resultp4,
1280 result_exponent,is_tiny);
1281 }
1282 Dbl_set_sign(resultp1,/*using*/sign_save);
1283 if (Dblext_isnotzero_mantissap3(resultp3) ||
1284 Dblext_isnotzero_mantissap4(resultp4)) {
1285 inexact = TRUE;
1286 switch(Rounding_mode()) {
1287 case ROUNDNEAREST: /* The default. */
1288 if (Dblext_isone_highp3(resultp3)) {
1289 /* at least 1/2 ulp */
1290 if (Dblext_isnotzero_low31p3(resultp3) ||
1291 Dblext_isnotzero_mantissap4(resultp4) ||
1292 Dblext_isone_lowp2(resultp2)) {
1293 /* either exactly half way and odd or
1294 * more than 1/2ulp */
1295 Dbl_increment(resultp1,resultp2);
1296 }
1297 }
1298 break;
1299
1300 case ROUNDPLUS:
1301 if (Dbl_iszero_sign(resultp1)) {
1302 /* Round up positive results */
1303 Dbl_increment(resultp1,resultp2);
1304 }
1305 break;
1306
1307 case ROUNDMINUS:
1308 if (Dbl_isone_sign(resultp1)) {
1309 /* Round down negative results */
1310 Dbl_increment(resultp1,resultp2);
1311 }
1312
1313 case ROUNDZERO:;
1314 /* truncate is simple */
1315 } /* end switch... */
1316 if (Dbl_isone_hiddenoverflow(resultp1)) result_exponent++;
1317 }
1318 if (result_exponent >= DBL_INFINITY_EXPONENT) {
1319 /* Overflow */
1320 if (Is_overflowtrap_enabled()) {
1321 /*
1322 * Adjust bias of result
1323 */
1324 Dbl_setwrapped_exponent(resultp1,result_exponent,ovfl);
1325 Dbl_copytoptr(resultp1,resultp2,dstptr);
1326 if (inexact)
1327 if (Is_inexacttrap_enabled())
1328 return (OPC_2E_OVERFLOWEXCEPTION |
1329 OPC_2E_INEXACTEXCEPTION);
1330 else Set_inexactflag();
1331 return (OPC_2E_OVERFLOWEXCEPTION);
1332 }
1333 inexact = TRUE;
1334 Set_overflowflag();
1335 Dbl_setoverflow(resultp1,resultp2);
1336 } else if (result_exponent <= 0) { /* underflow case */
1337 if (Is_underflowtrap_enabled()) {
1338 /*
1339 * Adjust bias of result
1340 */
1341 Dbl_setwrapped_exponent(resultp1,result_exponent,unfl);
1342 Dbl_copytoptr(resultp1,resultp2,dstptr);
1343 if (inexact)
1344 if (Is_inexacttrap_enabled())
1345 return (OPC_2E_UNDERFLOWEXCEPTION |
1346 OPC_2E_INEXACTEXCEPTION);
1347 else Set_inexactflag();
1348 return(OPC_2E_UNDERFLOWEXCEPTION);
1349 }
1350 else if (inexact && is_tiny) Set_underflowflag();
1351 }
1352 else Dbl_set_exponent(resultp1,result_exponent);
1353 Dbl_copytoptr(resultp1,resultp2,dstptr);
1354 if (inexact)
1355 if (Is_inexacttrap_enabled()) return(OPC_2E_INEXACTEXCEPTION);
1356 else Set_inexactflag();
1357 return(NOEXCEPTION);
1358 }
1359
1360 /*
1361 * Single Floating-point Multiply Fused Add
1362 */
1363
1364 sgl_fmpyfadd(src1ptr,src2ptr,src3ptr,status,dstptr)
1365
1366 sgl_floating_point *src1ptr, *src2ptr, *src3ptr, *dstptr;
1367 unsigned int *status;
1368 {
1369 unsigned int opnd1, opnd2, opnd3;
1370 register unsigned int tmpresp1, tmpresp2;
1371 unsigned int rightp1, rightp2;
1372 unsigned int resultp1, resultp2 = 0;
1373 register int mpy_exponent, add_exponent, count;
1374 boolean inexact = FALSE, is_tiny = FALSE;
1375
1376 unsigned int signlessleft1, signlessright1, save;
1377 register int result_exponent, diff_exponent;
1378 int sign_save, jumpsize;
1379
1380 Sgl_copyfromptr(src1ptr,opnd1);
1381 Sgl_copyfromptr(src2ptr,opnd2);
1382 Sgl_copyfromptr(src3ptr,opnd3);
1383
1384 /*
1385 * set sign bit of result of multiply
1386 */
1387 if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2))
1388 Sgl_setnegativezero(resultp1);
1389 else Sgl_setzero(resultp1);
1390
1391 /*
1392 * Generate multiply exponent
1393 */
1394 mpy_exponent = Sgl_exponent(opnd1) + Sgl_exponent(opnd2) - SGL_BIAS;
1395
1396 /*
1397 * check first operand for NaN's or infinity
1398 */
1399 if (Sgl_isinfinity_exponent(opnd1)) {
1400 if (Sgl_iszero_mantissa(opnd1)) {
1401 if (Sgl_isnotnan(opnd2) && Sgl_isnotnan(opnd3)) {
1402 if (Sgl_iszero_exponentmantissa(opnd2)) {
1403 /*
1404 * invalid since operands are infinity
1405 * and zero
1406 */
1407 if (Is_invalidtrap_enabled())
1408 return(OPC_2E_INVALIDEXCEPTION);
1409 Set_invalidflag();
1410 Sgl_makequietnan(resultp1);
1411 Sgl_copytoptr(resultp1,dstptr);
1412 return(NOEXCEPTION);
1413 }
1414 /*
1415 * Check third operand for infinity with a
1416 * sign opposite of the multiply result
1417 */
1418 if (Sgl_isinfinity(opnd3) &&
1419 (Sgl_sign(resultp1) ^ Sgl_sign(opnd3))) {
1420 /*
1421 * invalid since attempting a magnitude
1422 * subtraction of infinities
1423 */
1424 if (Is_invalidtrap_enabled())
1425 return(OPC_2E_INVALIDEXCEPTION);
1426 Set_invalidflag();
1427 Sgl_makequietnan(resultp1);
1428 Sgl_copytoptr(resultp1,dstptr);
1429 return(NOEXCEPTION);
1430 }
1431
1432 /*
1433 * return infinity
1434 */
1435 Sgl_setinfinity_exponentmantissa(resultp1);
1436 Sgl_copytoptr(resultp1,dstptr);
1437 return(NOEXCEPTION);
1438 }
1439 }
1440 else {
1441 /*
1442 * is NaN; signaling or quiet?
1443 */
1444 if (Sgl_isone_signaling(opnd1)) {
1445 /* trap if INVALIDTRAP enabled */
1446 if (Is_invalidtrap_enabled())
1447 return(OPC_2E_INVALIDEXCEPTION);
1448 /* make NaN quiet */
1449 Set_invalidflag();
1450 Sgl_set_quiet(opnd1);
1451 }
1452 /*
1453 * is second operand a signaling NaN?
1454 */
1455 else if (Sgl_is_signalingnan(opnd2)) {
1456 /* trap if INVALIDTRAP enabled */
1457 if (Is_invalidtrap_enabled())
1458 return(OPC_2E_INVALIDEXCEPTION);
1459 /* make NaN quiet */
1460 Set_invalidflag();
1461 Sgl_set_quiet(opnd2);
1462 Sgl_copytoptr(opnd2,dstptr);
1463 return(NOEXCEPTION);
1464 }
1465 /*
1466 * is third operand a signaling NaN?
1467 */
1468 else if (Sgl_is_signalingnan(opnd3)) {
1469 /* trap if INVALIDTRAP enabled */
1470 if (Is_invalidtrap_enabled())
1471 return(OPC_2E_INVALIDEXCEPTION);
1472 /* make NaN quiet */
1473 Set_invalidflag();
1474 Sgl_set_quiet(opnd3);
1475 Sgl_copytoptr(opnd3,dstptr);
1476 return(NOEXCEPTION);
1477 }
1478 /*
1479 * return quiet NaN
1480 */
1481 Sgl_copytoptr(opnd1,dstptr);
1482 return(NOEXCEPTION);
1483 }
1484 }
1485
1486 /*
1487 * check second operand for NaN's or infinity
1488 */
1489 if (Sgl_isinfinity_exponent(opnd2)) {
1490 if (Sgl_iszero_mantissa(opnd2)) {
1491 if (Sgl_isnotnan(opnd3)) {
1492 if (Sgl_iszero_exponentmantissa(opnd1)) {
1493 /*
1494 * invalid since multiply operands are
1495 * zero & infinity
1496 */
1497 if (Is_invalidtrap_enabled())
1498 return(OPC_2E_INVALIDEXCEPTION);
1499 Set_invalidflag();
1500 Sgl_makequietnan(opnd2);
1501 Sgl_copytoptr(opnd2,dstptr);
1502 return(NOEXCEPTION);
1503 }
1504
1505 /*
1506 * Check third operand for infinity with a
1507 * sign opposite of the multiply result
1508 */
1509 if (Sgl_isinfinity(opnd3) &&
1510 (Sgl_sign(resultp1) ^ Sgl_sign(opnd3))) {
1511 /*
1512 * invalid since attempting a magnitude
1513 * subtraction of infinities
1514 */
1515 if (Is_invalidtrap_enabled())
1516 return(OPC_2E_INVALIDEXCEPTION);
1517 Set_invalidflag();
1518 Sgl_makequietnan(resultp1);
1519 Sgl_copytoptr(resultp1,dstptr);
1520 return(NOEXCEPTION);
1521 }
1522
1523 /*
1524 * return infinity
1525 */
1526 Sgl_setinfinity_exponentmantissa(resultp1);
1527 Sgl_copytoptr(resultp1,dstptr);
1528 return(NOEXCEPTION);
1529 }
1530 }
1531 else {
1532 /*
1533 * is NaN; signaling or quiet?
1534 */
1535 if (Sgl_isone_signaling(opnd2)) {
1536 /* trap if INVALIDTRAP enabled */
1537 if (Is_invalidtrap_enabled())
1538 return(OPC_2E_INVALIDEXCEPTION);
1539 /* make NaN quiet */
1540 Set_invalidflag();
1541 Sgl_set_quiet(opnd2);
1542 }
1543 /*
1544 * is third operand a signaling NaN?
1545 */
1546 else if (Sgl_is_signalingnan(opnd3)) {
1547 /* trap if INVALIDTRAP enabled */
1548 if (Is_invalidtrap_enabled())
1549 return(OPC_2E_INVALIDEXCEPTION);
1550 /* make NaN quiet */
1551 Set_invalidflag();
1552 Sgl_set_quiet(opnd3);
1553 Sgl_copytoptr(opnd3,dstptr);
1554 return(NOEXCEPTION);
1555 }
1556 /*
1557 * return quiet NaN
1558 */
1559 Sgl_copytoptr(opnd2,dstptr);
1560 return(NOEXCEPTION);
1561 }
1562 }
1563
1564 /*
1565 * check third operand for NaN's or infinity
1566 */
1567 if (Sgl_isinfinity_exponent(opnd3)) {
1568 if (Sgl_iszero_mantissa(opnd3)) {
1569 /* return infinity */
1570 Sgl_copytoptr(opnd3,dstptr);
1571 return(NOEXCEPTION);
1572 } else {
1573 /*
1574 * is NaN; signaling or quiet?
1575 */
1576 if (Sgl_isone_signaling(opnd3)) {
1577 /* trap if INVALIDTRAP enabled */
1578 if (Is_invalidtrap_enabled())
1579 return(OPC_2E_INVALIDEXCEPTION);
1580 /* make NaN quiet */
1581 Set_invalidflag();
1582 Sgl_set_quiet(opnd3);
1583 }
1584 /*
1585 * return quiet NaN
1586 */
1587 Sgl_copytoptr(opnd3,dstptr);
1588 return(NOEXCEPTION);
1589 }
1590 }
1591
1592 /*
1593 * Generate multiply mantissa
1594 */
1595 if (Sgl_isnotzero_exponent(opnd1)) {
1596 /* set hidden bit */
1597 Sgl_clear_signexponent_set_hidden(opnd1);
1598 }
1599 else {
1600 /* check for zero */
1601 if (Sgl_iszero_mantissa(opnd1)) {
1602 /*
1603 * Perform the add opnd3 with zero here.
1604 */
1605 if (Sgl_iszero_exponentmantissa(opnd3)) {
1606 if (Is_rounding_mode(ROUNDMINUS)) {
1607 Sgl_or_signs(opnd3,resultp1);
1608 } else {
1609 Sgl_and_signs(opnd3,resultp1);
1610 }
1611 }
1612 /*
1613 * Now let's check for trapped underflow case.
1614 */
1615 else if (Sgl_iszero_exponent(opnd3) &&
1616 Is_underflowtrap_enabled()) {
1617 /* need to normalize results mantissa */
1618 sign_save = Sgl_signextendedsign(opnd3);
1619 result_exponent = 0;
1620 Sgl_leftshiftby1(opnd3);
1621 Sgl_normalize(opnd3,result_exponent);
1622 Sgl_set_sign(opnd3,/*using*/sign_save);
1623 Sgl_setwrapped_exponent(opnd3,result_exponent,
1624 unfl);
1625 Sgl_copytoptr(opnd3,dstptr);
1626 /* inexact = FALSE */
1627 return(OPC_2E_UNDERFLOWEXCEPTION);
1628 }
1629 Sgl_copytoptr(opnd3,dstptr);
1630 return(NOEXCEPTION);
1631 }
1632 /* is denormalized, adjust exponent */
1633 Sgl_clear_signexponent(opnd1);
1634 Sgl_leftshiftby1(opnd1);
1635 Sgl_normalize(opnd1,mpy_exponent);
1636 }
1637 /* opnd2 needs to have hidden bit set with msb in hidden bit */
1638 if (Sgl_isnotzero_exponent(opnd2)) {
1639 Sgl_clear_signexponent_set_hidden(opnd2);
1640 }
1641 else {
1642 /* check for zero */
1643 if (Sgl_iszero_mantissa(opnd2)) {
1644 /*
1645 * Perform the add opnd3 with zero here.
1646 */
1647 if (Sgl_iszero_exponentmantissa(opnd3)) {
1648 if (Is_rounding_mode(ROUNDMINUS)) {
1649 Sgl_or_signs(opnd3,resultp1);
1650 } else {
1651 Sgl_and_signs(opnd3,resultp1);
1652 }
1653 }
1654 /*
1655 * Now let's check for trapped underflow case.
1656 */
1657 else if (Sgl_iszero_exponent(opnd3) &&
1658 Is_underflowtrap_enabled()) {
1659 /* need to normalize results mantissa */
1660 sign_save = Sgl_signextendedsign(opnd3);
1661 result_exponent = 0;
1662 Sgl_leftshiftby1(opnd3);
1663 Sgl_normalize(opnd3,result_exponent);
1664 Sgl_set_sign(opnd3,/*using*/sign_save);
1665 Sgl_setwrapped_exponent(opnd3,result_exponent,
1666 unfl);
1667 Sgl_copytoptr(opnd3,dstptr);
1668 /* inexact = FALSE */
1669 return(OPC_2E_UNDERFLOWEXCEPTION);
1670 }
1671 Sgl_copytoptr(opnd3,dstptr);
1672 return(NOEXCEPTION);
1673 }
1674 /* is denormalized; want to normalize */
1675 Sgl_clear_signexponent(opnd2);
1676 Sgl_leftshiftby1(opnd2);
1677 Sgl_normalize(opnd2,mpy_exponent);
1678 }
1679
1680 /* Multiply the first two source mantissas together */
1681
1682 /*
1683 * The intermediate result will be kept in tmpres,
1684 * which needs enough room for 106 bits of mantissa,
1685 * so lets call it a Double extended.
1686 */
1687 Sglext_setzero(tmpresp1,tmpresp2);
1688
1689 /*
1690 * Four bits at a time are inspected in each loop, and a
1691 * simple shift and add multiply algorithm is used.
1692 */
1693 for (count = SGL_P-1; count >= 0; count -= 4) {
1694 Sglext_rightshiftby4(tmpresp1,tmpresp2);
1695 if (Sbit28(opnd1)) {
1696 /* Twoword_add should be an ADD followed by 2 ADDC's */
1697 Twoword_add(tmpresp1, tmpresp2, opnd2<<3, 0);
1698 }
1699 if (Sbit29(opnd1)) {
1700 Twoword_add(tmpresp1, tmpresp2, opnd2<<2, 0);
1701 }
1702 if (Sbit30(opnd1)) {
1703 Twoword_add(tmpresp1, tmpresp2, opnd2<<1, 0);
1704 }
1705 if (Sbit31(opnd1)) {
1706 Twoword_add(tmpresp1, tmpresp2, opnd2, 0);
1707 }
1708 Sgl_rightshiftby4(opnd1);
1709 }
1710 if (Is_sexthiddenoverflow(tmpresp1)) {
1711 /* result mantissa >= 2 (mantissa overflow) */
1712 mpy_exponent++;
1713 Sglext_rightshiftby4(tmpresp1,tmpresp2);
1714 } else {
1715 Sglext_rightshiftby3(tmpresp1,tmpresp2);
1716 }
1717
1718 /*
1719 * Restore the sign of the mpy result which was saved in resultp1.
1720 * The exponent will continue to be kept in mpy_exponent.
1721 */
1722 Sglext_set_sign(tmpresp1,Sgl_sign(resultp1));
1723
1724 /*
1725 * No rounding is required, since the result of the multiply
1726 * is exact in the extended format.
1727 */
1728
1729 /*
1730 * Now we are ready to perform the add portion of the operation.
1731 *
1732 * The exponents need to be kept as integers for now, since the
1733 * multiply result might not fit into the exponent field. We
1734 * can't overflow or underflow because of this yet, since the
1735 * add could bring the final result back into range.
1736 */
1737 add_exponent = Sgl_exponent(opnd3);
1738
1739 /*
1740 * Check for denormalized or zero add operand.
1741 */
1742 if (add_exponent == 0) {
1743 /* check for zero */
1744 if (Sgl_iszero_mantissa(opnd3)) {
1745 /* right is zero */
1746 /* Left can't be zero and must be result.
1747 *
1748 * The final result is now in tmpres and mpy_exponent,
1749 * and needs to be rounded and squeezed back into
1750 * double precision format from double extended.
1751 */
1752 result_exponent = mpy_exponent;
1753 Sglext_copy(tmpresp1,tmpresp2,resultp1,resultp2);
1754 sign_save = Sgl_signextendedsign(resultp1);/*save sign*/
1755 goto round;
1756 }
1757
1758 /*
1759 * Neither are zeroes.
1760 * Adjust exponent and normalize add operand.
1761 */
1762 sign_save = Sgl_signextendedsign(opnd3); /* save sign */
1763 Sgl_clear_signexponent(opnd3);
1764 Sgl_leftshiftby1(opnd3);
1765 Sgl_normalize(opnd3,add_exponent);
1766 Sgl_set_sign(opnd3,sign_save); /* restore sign */
1767 } else {
1768 Sgl_clear_exponent_set_hidden(opnd3);
1769 }
1770 /*
1771 * Copy opnd3 to the double extended variable called right.
1772 */
1773 Sgl_copyto_sglext(opnd3,rightp1,rightp2);
1774
1775 /*
1776 * A zero "save" helps discover equal operands (for later),
1777 * and is used in swapping operands (if needed).
1778 */
1779 Sglext_xortointp1(tmpresp1,rightp1,/*to*/save);
1780
1781 /*
1782 * Compare magnitude of operands.
1783 */
1784 Sglext_copytoint_exponentmantissa(tmpresp1,signlessleft1);
1785 Sglext_copytoint_exponentmantissa(rightp1,signlessright1);
1786 if (mpy_exponent < add_exponent || mpy_exponent == add_exponent &&
1787 Sglext_ismagnitudeless(signlessleft1,signlessright1)) {
1788 /*
1789 * Set the left operand to the larger one by XOR swap.
1790 * First finish the first word "save".
1791 */
1792 Sglext_xorfromintp1(save,rightp1,/*to*/rightp1);
1793 Sglext_xorfromintp1(save,tmpresp1,/*to*/tmpresp1);
1794 Sglext_swap_lower(tmpresp2,rightp2);
1795 /* also setup exponents used in rest of routine */
1796 diff_exponent = add_exponent - mpy_exponent;
1797 result_exponent = add_exponent;
1798 } else {
1799 /* also setup exponents used in rest of routine */
1800 diff_exponent = mpy_exponent - add_exponent;
1801 result_exponent = mpy_exponent;
1802 }
1803 /* Invariant: left is not smaller than right. */
1804
1805 /*
1806 * Special case alignment of operands that would force alignment
1807 * beyond the extent of the extension. A further optimization
1808 * could special case this but only reduces the path length for
1809 * this infrequent case.
1810 */
1811 if (diff_exponent > SGLEXT_THRESHOLD) {
1812 diff_exponent = SGLEXT_THRESHOLD;
1813 }
1814
1815 /* Align right operand by shifting it to the right */
1816 Sglext_clear_sign(rightp1);
1817 Sglext_right_align(rightp1,rightp2,/*shifted by*/diff_exponent);
1818
1819 /* Treat sum and difference of the operands separately. */
1820 if ((int)save < 0) {
1821 /*
1822 * Difference of the two operands. Overflow can occur if the
1823 * multiply overflowed. A borrow can occur out of the hidden
1824 * bit and force a post normalization phase.
1825 */
1826 Sglext_subtract(tmpresp1,tmpresp2, rightp1,rightp2,
1827 resultp1,resultp2);
1828 sign_save = Sgl_signextendedsign(resultp1);
1829 if (Sgl_iszero_hidden(resultp1)) {
1830 /* Handle normalization */
1831 /* A straightforward algorithm would now shift the
1832 * result and extension left until the hidden bit
1833 * becomes one. Not all of the extension bits need
1834 * participate in the shift. Only the two most
1835 * significant bits (round and guard) are needed.
1836 * If only a single shift is needed then the guard
1837 * bit becomes a significant low order bit and the
1838 * extension must participate in the rounding.
1839 * If more than a single shift is needed, then all
1840 * bits to the right of the guard bit are zeros,
1841 * and the guard bit may or may not be zero. */
1842 Sglext_leftshiftby1(resultp1,resultp2);
1843
1844 /* Need to check for a zero result. The sign and
1845 * exponent fields have already been zeroed. The more
1846 * efficient test of the full object can be used.
1847 */
1848 if (Sglext_iszero(resultp1,resultp2)) {
1849 /* Must have been "x-x" or "x+(-x)". */
1850 if (Is_rounding_mode(ROUNDMINUS))
1851 Sgl_setone_sign(resultp1);
1852 Sgl_copytoptr(resultp1,dstptr);
1853 return(NOEXCEPTION);
1854 }
1855 result_exponent--;
1856
1857 /* Look to see if normalization is finished. */
1858 if (Sgl_isone_hidden(resultp1)) {
1859 /* No further normalization is needed */
1860 goto round;
1861 }
1862
1863 /* Discover first one bit to determine shift amount.
1864 * Use a modified binary search. We have already
1865 * shifted the result one position right and still
1866 * not found a one so the remainder of the extension
1867 * must be zero and simplifies rounding. */
1868 /* Scan bytes */
1869 while (Sgl_iszero_hiddenhigh7mantissa(resultp1)) {
1870 Sglext_leftshiftby8(resultp1,resultp2);
1871 result_exponent -= 8;
1872 }
1873 /* Now narrow it down to the nibble */
1874 if (Sgl_iszero_hiddenhigh3mantissa(resultp1)) {
1875 /* The lower nibble contains the
1876 * normalizing one */
1877 Sglext_leftshiftby4(resultp1,resultp2);
1878 result_exponent -= 4;
1879 }
1880 /* Select case where first bit is set (already
1881 * normalized) otherwise select the proper shift. */
1882 jumpsize = Sgl_hiddenhigh3mantissa(resultp1);
1883 if (jumpsize <= 7) switch(jumpsize) {
1884 case 1:
1885 Sglext_leftshiftby3(resultp1,resultp2);
1886 result_exponent -= 3;
1887 break;
1888 case 2:
1889 case 3:
1890 Sglext_leftshiftby2(resultp1,resultp2);
1891 result_exponent -= 2;
1892 break;
1893 case 4:
1894 case 5:
1895 case 6:
1896 case 7:
1897 Sglext_leftshiftby1(resultp1,resultp2);
1898 result_exponent -= 1;
1899 break;
1900 }
1901 } /* end if (hidden...)... */
1902 /* Fall through and round */
1903 } /* end if (save < 0)... */
1904 else {
1905 /* Add magnitudes */
1906 Sglext_addition(tmpresp1,tmpresp2,
1907 rightp1,rightp2, /*to*/resultp1,resultp2);
1908 sign_save = Sgl_signextendedsign(resultp1);
1909 if (Sgl_isone_hiddenoverflow(resultp1)) {
1910 /* Prenormalization required. */
1911 Sglext_arithrightshiftby1(resultp1,resultp2);
1912 result_exponent++;
1913 } /* end if hiddenoverflow... */
1914 } /* end else ...add magnitudes... */
1915
1916 /* Round the result. If the extension and lower two words are
1917 * all zeros, then the result is exact. Otherwise round in the
1918 * correct direction. Underflow is possible. If a postnormalization
1919 * is necessary, then the mantissa is all zeros so no shift is needed.
1920 */
1921 round:
1922 if (result_exponent <= 0 && !Is_underflowtrap_enabled()) {
1923 Sglext_denormalize(resultp1,resultp2,result_exponent,is_tiny);
1924 }
1925 Sgl_set_sign(resultp1,/*using*/sign_save);
1926 if (Sglext_isnotzero_mantissap2(resultp2)) {
1927 inexact = TRUE;
1928 switch(Rounding_mode()) {
1929 case ROUNDNEAREST: /* The default. */
1930 if (Sglext_isone_highp2(resultp2)) {
1931 /* at least 1/2 ulp */
1932 if (Sglext_isnotzero_low31p2(resultp2) ||
1933 Sglext_isone_lowp1(resultp1)) {
1934 /* either exactly half way and odd or
1935 * more than 1/2ulp */
1936 Sgl_increment(resultp1);
1937 }
1938 }
1939 break;
1940
1941 case ROUNDPLUS:
1942 if (Sgl_iszero_sign(resultp1)) {
1943 /* Round up positive results */
1944 Sgl_increment(resultp1);
1945 }
1946 break;
1947
1948 case ROUNDMINUS:
1949 if (Sgl_isone_sign(resultp1)) {
1950 /* Round down negative results */
1951 Sgl_increment(resultp1);
1952 }
1953
1954 case ROUNDZERO:;
1955 /* truncate is simple */
1956 } /* end switch... */
1957 if (Sgl_isone_hiddenoverflow(resultp1)) result_exponent++;
1958 }
1959 if (result_exponent >= SGL_INFINITY_EXPONENT) {
1960 /* Overflow */
1961 if (Is_overflowtrap_enabled()) {
1962 /*
1963 * Adjust bias of result
1964 */
1965 Sgl_setwrapped_exponent(resultp1,result_exponent,ovfl);
1966 Sgl_copytoptr(resultp1,dstptr);
1967 if (inexact)
1968 if (Is_inexacttrap_enabled())
1969 return (OPC_2E_OVERFLOWEXCEPTION |
1970 OPC_2E_INEXACTEXCEPTION);
1971 else Set_inexactflag();
1972 return (OPC_2E_OVERFLOWEXCEPTION);
1973 }
1974 inexact = TRUE;
1975 Set_overflowflag();
1976 Sgl_setoverflow(resultp1);
1977 } else if (result_exponent <= 0) { /* underflow case */
1978 if (Is_underflowtrap_enabled()) {
1979 /*
1980 * Adjust bias of result
1981 */
1982 Sgl_setwrapped_exponent(resultp1,result_exponent,unfl);
1983 Sgl_copytoptr(resultp1,dstptr);
1984 if (inexact)
1985 if (Is_inexacttrap_enabled())
1986 return (OPC_2E_UNDERFLOWEXCEPTION |
1987 OPC_2E_INEXACTEXCEPTION);
1988 else Set_inexactflag();
1989 return(OPC_2E_UNDERFLOWEXCEPTION);
1990 }
1991 else if (inexact && is_tiny) Set_underflowflag();
1992 }
1993 else Sgl_set_exponent(resultp1,result_exponent);
1994 Sgl_copytoptr(resultp1,dstptr);
1995 if (inexact)
1996 if (Is_inexacttrap_enabled()) return(OPC_2E_INEXACTEXCEPTION);
1997 else Set_inexactflag();
1998 return(NOEXCEPTION);
1999 }
2000
2001 /*
2002 * Single Floating-point Multiply Negate Fused Add
2003 */
2004
2005 sgl_fmpynfadd(src1ptr,src2ptr,src3ptr,status,dstptr)
2006
2007 sgl_floating_point *src1ptr, *src2ptr, *src3ptr, *dstptr;
2008 unsigned int *status;
2009 {
2010 unsigned int opnd1, opnd2, opnd3;
2011 register unsigned int tmpresp1, tmpresp2;
2012 unsigned int rightp1, rightp2;
2013 unsigned int resultp1, resultp2 = 0;
2014 register int mpy_exponent, add_exponent, count;
2015 boolean inexact = FALSE, is_tiny = FALSE;
2016
2017 unsigned int signlessleft1, signlessright1, save;
2018 register int result_exponent, diff_exponent;
2019 int sign_save, jumpsize;
2020
2021 Sgl_copyfromptr(src1ptr,opnd1);
2022 Sgl_copyfromptr(src2ptr,opnd2);
2023 Sgl_copyfromptr(src3ptr,opnd3);
2024
2025 /*
2026 * set sign bit of result of multiply
2027 */
2028 if (Sgl_sign(opnd1) ^ Sgl_sign(opnd2))
2029 Sgl_setzero(resultp1);
2030 else
2031 Sgl_setnegativezero(resultp1);
2032
2033 /*
2034 * Generate multiply exponent
2035 */
2036 mpy_exponent = Sgl_exponent(opnd1) + Sgl_exponent(opnd2) - SGL_BIAS;
2037
2038 /*
2039 * check first operand for NaN's or infinity
2040 */
2041 if (Sgl_isinfinity_exponent(opnd1)) {
2042 if (Sgl_iszero_mantissa(opnd1)) {
2043 if (Sgl_isnotnan(opnd2) && Sgl_isnotnan(opnd3)) {
2044 if (Sgl_iszero_exponentmantissa(opnd2)) {
2045 /*
2046 * invalid since operands are infinity
2047 * and zero
2048 */
2049 if (Is_invalidtrap_enabled())
2050 return(OPC_2E_INVALIDEXCEPTION);
2051 Set_invalidflag();
2052 Sgl_makequietnan(resultp1);
2053 Sgl_copytoptr(resultp1,dstptr);
2054 return(NOEXCEPTION);
2055 }
2056 /*
2057 * Check third operand for infinity with a
2058 * sign opposite of the multiply result
2059 */
2060 if (Sgl_isinfinity(opnd3) &&
2061 (Sgl_sign(resultp1) ^ Sgl_sign(opnd3))) {
2062 /*
2063 * invalid since attempting a magnitude
2064 * subtraction of infinities
2065 */
2066 if (Is_invalidtrap_enabled())
2067 return(OPC_2E_INVALIDEXCEPTION);
2068 Set_invalidflag();
2069 Sgl_makequietnan(resultp1);
2070 Sgl_copytoptr(resultp1,dstptr);
2071 return(NOEXCEPTION);
2072 }
2073
2074 /*
2075 * return infinity
2076 */
2077 Sgl_setinfinity_exponentmantissa(resultp1);
2078 Sgl_copytoptr(resultp1,dstptr);
2079 return(NOEXCEPTION);
2080 }
2081 }
2082 else {
2083 /*
2084 * is NaN; signaling or quiet?
2085 */
2086 if (Sgl_isone_signaling(opnd1)) {
2087 /* trap if INVALIDTRAP enabled */
2088 if (Is_invalidtrap_enabled())
2089 return(OPC_2E_INVALIDEXCEPTION);
2090 /* make NaN quiet */
2091 Set_invalidflag();
2092 Sgl_set_quiet(opnd1);
2093 }
2094 /*
2095 * is second operand a signaling NaN?
2096 */
2097 else if (Sgl_is_signalingnan(opnd2)) {
2098 /* trap if INVALIDTRAP enabled */
2099 if (Is_invalidtrap_enabled())
2100 return(OPC_2E_INVALIDEXCEPTION);
2101 /* make NaN quiet */
2102 Set_invalidflag();
2103 Sgl_set_quiet(opnd2);
2104 Sgl_copytoptr(opnd2,dstptr);
2105 return(NOEXCEPTION);
2106 }
2107 /*
2108 * is third operand a signaling NaN?
2109 */
2110 else if (Sgl_is_signalingnan(opnd3)) {
2111 /* trap if INVALIDTRAP enabled */
2112 if (Is_invalidtrap_enabled())
2113 return(OPC_2E_INVALIDEXCEPTION);
2114 /* make NaN quiet */
2115 Set_invalidflag();
2116 Sgl_set_quiet(opnd3);
2117 Sgl_copytoptr(opnd3,dstptr);
2118 return(NOEXCEPTION);
2119 }
2120 /*
2121 * return quiet NaN
2122 */
2123 Sgl_copytoptr(opnd1,dstptr);
2124 return(NOEXCEPTION);
2125 }
2126 }
2127
2128 /*
2129 * check second operand for NaN's or infinity
2130 */
2131 if (Sgl_isinfinity_exponent(opnd2)) {
2132 if (Sgl_iszero_mantissa(opnd2)) {
2133 if (Sgl_isnotnan(opnd3)) {
2134 if (Sgl_iszero_exponentmantissa(opnd1)) {
2135 /*
2136 * invalid since multiply operands are
2137 * zero & infinity
2138 */
2139 if (Is_invalidtrap_enabled())
2140 return(OPC_2E_INVALIDEXCEPTION);
2141 Set_invalidflag();
2142 Sgl_makequietnan(opnd2);
2143 Sgl_copytoptr(opnd2,dstptr);
2144 return(NOEXCEPTION);
2145 }
2146
2147 /*
2148 * Check third operand for infinity with a
2149 * sign opposite of the multiply result
2150 */
2151 if (Sgl_isinfinity(opnd3) &&
2152 (Sgl_sign(resultp1) ^ Sgl_sign(opnd3))) {
2153 /*
2154 * invalid since attempting a magnitude
2155 * subtraction of infinities
2156 */
2157 if (Is_invalidtrap_enabled())
2158 return(OPC_2E_INVALIDEXCEPTION);
2159 Set_invalidflag();
2160 Sgl_makequietnan(resultp1);
2161 Sgl_copytoptr(resultp1,dstptr);
2162 return(NOEXCEPTION);
2163 }
2164
2165 /*
2166 * return infinity
2167 */
2168 Sgl_setinfinity_exponentmantissa(resultp1);
2169 Sgl_copytoptr(resultp1,dstptr);
2170 return(NOEXCEPTION);
2171 }
2172 }
2173 else {
2174 /*
2175 * is NaN; signaling or quiet?
2176 */
2177 if (Sgl_isone_signaling(opnd2)) {
2178 /* trap if INVALIDTRAP enabled */
2179 if (Is_invalidtrap_enabled())
2180 return(OPC_2E_INVALIDEXCEPTION);
2181 /* make NaN quiet */
2182 Set_invalidflag();
2183 Sgl_set_quiet(opnd2);
2184 }
2185 /*
2186 * is third operand a signaling NaN?
2187 */
2188 else if (Sgl_is_signalingnan(opnd3)) {
2189 /* trap if INVALIDTRAP enabled */
2190 if (Is_invalidtrap_enabled())
2191 return(OPC_2E_INVALIDEXCEPTION);
2192 /* make NaN quiet */
2193 Set_invalidflag();
2194 Sgl_set_quiet(opnd3);
2195 Sgl_copytoptr(opnd3,dstptr);
2196 return(NOEXCEPTION);
2197 }
2198 /*
2199 * return quiet NaN
2200 */
2201 Sgl_copytoptr(opnd2,dstptr);
2202 return(NOEXCEPTION);
2203 }
2204 }
2205
2206 /*
2207 * check third operand for NaN's or infinity
2208 */
2209 if (Sgl_isinfinity_exponent(opnd3)) {
2210 if (Sgl_iszero_mantissa(opnd3)) {
2211 /* return infinity */
2212 Sgl_copytoptr(opnd3,dstptr);
2213 return(NOEXCEPTION);
2214 } else {
2215 /*
2216 * is NaN; signaling or quiet?
2217 */
2218 if (Sgl_isone_signaling(opnd3)) {
2219 /* trap if INVALIDTRAP enabled */
2220 if (Is_invalidtrap_enabled())
2221 return(OPC_2E_INVALIDEXCEPTION);
2222 /* make NaN quiet */
2223 Set_invalidflag();
2224 Sgl_set_quiet(opnd3);
2225 }
2226 /*
2227 * return quiet NaN
2228 */
2229 Sgl_copytoptr(opnd3,dstptr);
2230 return(NOEXCEPTION);
2231 }
2232 }
2233
2234 /*
2235 * Generate multiply mantissa
2236 */
2237 if (Sgl_isnotzero_exponent(opnd1)) {
2238 /* set hidden bit */
2239 Sgl_clear_signexponent_set_hidden(opnd1);
2240 }
2241 else {
2242 /* check for zero */
2243 if (Sgl_iszero_mantissa(opnd1)) {
2244 /*
2245 * Perform the add opnd3 with zero here.
2246 */
2247 if (Sgl_iszero_exponentmantissa(opnd3)) {
2248 if (Is_rounding_mode(ROUNDMINUS)) {
2249 Sgl_or_signs(opnd3,resultp1);
2250 } else {
2251 Sgl_and_signs(opnd3,resultp1);
2252 }
2253 }
2254 /*
2255 * Now let's check for trapped underflow case.
2256 */
2257 else if (Sgl_iszero_exponent(opnd3) &&
2258 Is_underflowtrap_enabled()) {
2259 /* need to normalize results mantissa */
2260 sign_save = Sgl_signextendedsign(opnd3);
2261 result_exponent = 0;
2262 Sgl_leftshiftby1(opnd3);
2263 Sgl_normalize(opnd3,result_exponent);
2264 Sgl_set_sign(opnd3,/*using*/sign_save);
2265 Sgl_setwrapped_exponent(opnd3,result_exponent,
2266 unfl);
2267 Sgl_copytoptr(opnd3,dstptr);
2268 /* inexact = FALSE */
2269 return(OPC_2E_UNDERFLOWEXCEPTION);
2270 }
2271 Sgl_copytoptr(opnd3,dstptr);
2272 return(NOEXCEPTION);
2273 }
2274 /* is denormalized, adjust exponent */
2275 Sgl_clear_signexponent(opnd1);
2276 Sgl_leftshiftby1(opnd1);
2277 Sgl_normalize(opnd1,mpy_exponent);
2278 }
2279 /* opnd2 needs to have hidden bit set with msb in hidden bit */
2280 if (Sgl_isnotzero_exponent(opnd2)) {
2281 Sgl_clear_signexponent_set_hidden(opnd2);
2282 }
2283 else {
2284 /* check for zero */
2285 if (Sgl_iszero_mantissa(opnd2)) {
2286 /*
2287 * Perform the add opnd3 with zero here.
2288 */
2289 if (Sgl_iszero_exponentmantissa(opnd3)) {
2290 if (Is_rounding_mode(ROUNDMINUS)) {
2291 Sgl_or_signs(opnd3,resultp1);
2292 } else {
2293 Sgl_and_signs(opnd3,resultp1);
2294 }
2295 }
2296 /*
2297 * Now let's check for trapped underflow case.
2298 */
2299 else if (Sgl_iszero_exponent(opnd3) &&
2300 Is_underflowtrap_enabled()) {
2301 /* need to normalize results mantissa */
2302 sign_save = Sgl_signextendedsign(opnd3);
2303 result_exponent = 0;
2304 Sgl_leftshiftby1(opnd3);
2305 Sgl_normalize(opnd3,result_exponent);
2306 Sgl_set_sign(opnd3,/*using*/sign_save);
2307 Sgl_setwrapped_exponent(opnd3,result_exponent,
2308 unfl);
2309 Sgl_copytoptr(opnd3,dstptr);
2310 /* inexact = FALSE */
2311 return(OPC_2E_UNDERFLOWEXCEPTION);
2312 }
2313 Sgl_copytoptr(opnd3,dstptr);
2314 return(NOEXCEPTION);
2315 }
2316 /* is denormalized; want to normalize */
2317 Sgl_clear_signexponent(opnd2);
2318 Sgl_leftshiftby1(opnd2);
2319 Sgl_normalize(opnd2,mpy_exponent);
2320 }
2321
2322 /* Multiply the first two source mantissas together */
2323
2324 /*
2325 * The intermediate result will be kept in tmpres,
2326 * which needs enough room for 106 bits of mantissa,
2327 * so lets call it a Double extended.
2328 */
2329 Sglext_setzero(tmpresp1,tmpresp2);
2330
2331 /*
2332 * Four bits at a time are inspected in each loop, and a
2333 * simple shift and add multiply algorithm is used.
2334 */
2335 for (count = SGL_P-1; count >= 0; count -= 4) {
2336 Sglext_rightshiftby4(tmpresp1,tmpresp2);
2337 if (Sbit28(opnd1)) {
2338 /* Twoword_add should be an ADD followed by 2 ADDC's */
2339 Twoword_add(tmpresp1, tmpresp2, opnd2<<3, 0);
2340 }
2341 if (Sbit29(opnd1)) {
2342 Twoword_add(tmpresp1, tmpresp2, opnd2<<2, 0);
2343 }
2344 if (Sbit30(opnd1)) {
2345 Twoword_add(tmpresp1, tmpresp2, opnd2<<1, 0);
2346 }
2347 if (Sbit31(opnd1)) {
2348 Twoword_add(tmpresp1, tmpresp2, opnd2, 0);
2349 }
2350 Sgl_rightshiftby4(opnd1);
2351 }
2352 if (Is_sexthiddenoverflow(tmpresp1)) {
2353 /* result mantissa >= 2 (mantissa overflow) */
2354 mpy_exponent++;
2355 Sglext_rightshiftby4(tmpresp1,tmpresp2);
2356 } else {
2357 Sglext_rightshiftby3(tmpresp1,tmpresp2);
2358 }
2359
2360 /*
2361 * Restore the sign of the mpy result which was saved in resultp1.
2362 * The exponent will continue to be kept in mpy_exponent.
2363 */
2364 Sglext_set_sign(tmpresp1,Sgl_sign(resultp1));
2365
2366 /*
2367 * No rounding is required, since the result of the multiply
2368 * is exact in the extended format.
2369 */
2370
2371 /*
2372 * Now we are ready to perform the add portion of the operation.
2373 *
2374 * The exponents need to be kept as integers for now, since the
2375 * multiply result might not fit into the exponent field. We
2376 * can't overflow or underflow because of this yet, since the
2377 * add could bring the final result back into range.
2378 */
2379 add_exponent = Sgl_exponent(opnd3);
2380
2381 /*
2382 * Check for denormalized or zero add operand.
2383 */
2384 if (add_exponent == 0) {
2385 /* check for zero */
2386 if (Sgl_iszero_mantissa(opnd3)) {
2387 /* right is zero */
2388 /* Left can't be zero and must be result.
2389 *
2390 * The final result is now in tmpres and mpy_exponent,
2391 * and needs to be rounded and squeezed back into
2392 * double precision format from double extended.
2393 */
2394 result_exponent = mpy_exponent;
2395 Sglext_copy(tmpresp1,tmpresp2,resultp1,resultp2);
2396 sign_save = Sgl_signextendedsign(resultp1);/*save sign*/
2397 goto round;
2398 }
2399
2400 /*
2401 * Neither are zeroes.
2402 * Adjust exponent and normalize add operand.
2403 */
2404 sign_save = Sgl_signextendedsign(opnd3); /* save sign */
2405 Sgl_clear_signexponent(opnd3);
2406 Sgl_leftshiftby1(opnd3);
2407 Sgl_normalize(opnd3,add_exponent);
2408 Sgl_set_sign(opnd3,sign_save); /* restore sign */
2409 } else {
2410 Sgl_clear_exponent_set_hidden(opnd3);
2411 }
2412 /*
2413 * Copy opnd3 to the double extended variable called right.
2414 */
2415 Sgl_copyto_sglext(opnd3,rightp1,rightp2);
2416
2417 /*
2418 * A zero "save" helps discover equal operands (for later),
2419 * and is used in swapping operands (if needed).
2420 */
2421 Sglext_xortointp1(tmpresp1,rightp1,/*to*/save);
2422
2423 /*
2424 * Compare magnitude of operands.
2425 */
2426 Sglext_copytoint_exponentmantissa(tmpresp1,signlessleft1);
2427 Sglext_copytoint_exponentmantissa(rightp1,signlessright1);
2428 if (mpy_exponent < add_exponent || mpy_exponent == add_exponent &&
2429 Sglext_ismagnitudeless(signlessleft1,signlessright1)) {
2430 /*
2431 * Set the left operand to the larger one by XOR swap.
2432 * First finish the first word "save".
2433 */
2434 Sglext_xorfromintp1(save,rightp1,/*to*/rightp1);
2435 Sglext_xorfromintp1(save,tmpresp1,/*to*/tmpresp1);
2436 Sglext_swap_lower(tmpresp2,rightp2);
2437 /* also setup exponents used in rest of routine */
2438 diff_exponent = add_exponent - mpy_exponent;
2439 result_exponent = add_exponent;
2440 } else {
2441 /* also setup exponents used in rest of routine */
2442 diff_exponent = mpy_exponent - add_exponent;
2443 result_exponent = mpy_exponent;
2444 }
2445 /* Invariant: left is not smaller than right. */
2446
2447 /*
2448 * Special case alignment of operands that would force alignment
2449 * beyond the extent of the extension. A further optimization
2450 * could special case this but only reduces the path length for
2451 * this infrequent case.
2452 */
2453 if (diff_exponent > SGLEXT_THRESHOLD) {
2454 diff_exponent = SGLEXT_THRESHOLD;
2455 }
2456
2457 /* Align right operand by shifting it to the right */
2458 Sglext_clear_sign(rightp1);
2459 Sglext_right_align(rightp1,rightp2,/*shifted by*/diff_exponent);
2460
2461 /* Treat sum and difference of the operands separately. */
2462 if ((int)save < 0) {
2463 /*
2464 * Difference of the two operands. Overflow can occur if the
2465 * multiply overflowed. A borrow can occur out of the hidden
2466 * bit and force a post normalization phase.
2467 */
2468 Sglext_subtract(tmpresp1,tmpresp2, rightp1,rightp2,
2469 resultp1,resultp2);
2470 sign_save = Sgl_signextendedsign(resultp1);
2471 if (Sgl_iszero_hidden(resultp1)) {
2472 /* Handle normalization */
2473 /* A straightforward algorithm would now shift the
2474 * result and extension left until the hidden bit
2475 * becomes one. Not all of the extension bits need
2476 * participate in the shift. Only the two most
2477 * significant bits (round and guard) are needed.
2478 * If only a single shift is needed then the guard
2479 * bit becomes a significant low order bit and the
2480 * extension must participate in the rounding.
2481 * If more than a single shift is needed, then all
2482 * bits to the right of the guard bit are zeros,
2483 * and the guard bit may or may not be zero. */
2484 Sglext_leftshiftby1(resultp1,resultp2);
2485
2486 /* Need to check for a zero result. The sign and
2487 * exponent fields have already been zeroed. The more
2488 * efficient test of the full object can be used.
2489 */
2490 if (Sglext_iszero(resultp1,resultp2)) {
2491 /* Must have been "x-x" or "x+(-x)". */
2492 if (Is_rounding_mode(ROUNDMINUS))
2493 Sgl_setone_sign(resultp1);
2494 Sgl_copytoptr(resultp1,dstptr);
2495 return(NOEXCEPTION);
2496 }
2497 result_exponent--;
2498
2499 /* Look to see if normalization is finished. */
2500 if (Sgl_isone_hidden(resultp1)) {
2501 /* No further normalization is needed */
2502 goto round;
2503 }
2504
2505 /* Discover first one bit to determine shift amount.
2506 * Use a modified binary search. We have already
2507 * shifted the result one position right and still
2508 * not found a one so the remainder of the extension
2509 * must be zero and simplifies rounding. */
2510 /* Scan bytes */
2511 while (Sgl_iszero_hiddenhigh7mantissa(resultp1)) {
2512 Sglext_leftshiftby8(resultp1,resultp2);
2513 result_exponent -= 8;
2514 }
2515 /* Now narrow it down to the nibble */
2516 if (Sgl_iszero_hiddenhigh3mantissa(resultp1)) {
2517 /* The lower nibble contains the
2518 * normalizing one */
2519 Sglext_leftshiftby4(resultp1,resultp2);
2520 result_exponent -= 4;
2521 }
2522 /* Select case where first bit is set (already
2523 * normalized) otherwise select the proper shift. */
2524 jumpsize = Sgl_hiddenhigh3mantissa(resultp1);
2525 if (jumpsize <= 7) switch(jumpsize) {
2526 case 1:
2527 Sglext_leftshiftby3(resultp1,resultp2);
2528 result_exponent -= 3;
2529 break;
2530 case 2:
2531 case 3:
2532 Sglext_leftshiftby2(resultp1,resultp2);
2533 result_exponent -= 2;
2534 break;
2535 case 4:
2536 case 5:
2537 case 6:
2538 case 7:
2539 Sglext_leftshiftby1(resultp1,resultp2);
2540 result_exponent -= 1;
2541 break;
2542 }
2543 } /* end if (hidden...)... */
2544 /* Fall through and round */
2545 } /* end if (save < 0)... */
2546 else {
2547 /* Add magnitudes */
2548 Sglext_addition(tmpresp1,tmpresp2,
2549 rightp1,rightp2, /*to*/resultp1,resultp2);
2550 sign_save = Sgl_signextendedsign(resultp1);
2551 if (Sgl_isone_hiddenoverflow(resultp1)) {
2552 /* Prenormalization required. */
2553 Sglext_arithrightshiftby1(resultp1,resultp2);
2554 result_exponent++;
2555 } /* end if hiddenoverflow... */
2556 } /* end else ...add magnitudes... */
2557
2558 /* Round the result. If the extension and lower two words are
2559 * all zeros, then the result is exact. Otherwise round in the
2560 * correct direction. Underflow is possible. If a postnormalization
2561 * is necessary, then the mantissa is all zeros so no shift is needed.
2562 */
2563 round:
2564 if (result_exponent <= 0 && !Is_underflowtrap_enabled()) {
2565 Sglext_denormalize(resultp1,resultp2,result_exponent,is_tiny);
2566 }
2567 Sgl_set_sign(resultp1,/*using*/sign_save);
2568 if (Sglext_isnotzero_mantissap2(resultp2)) {
2569 inexact = TRUE;
2570 switch(Rounding_mode()) {
2571 case ROUNDNEAREST: /* The default. */
2572 if (Sglext_isone_highp2(resultp2)) {
2573 /* at least 1/2 ulp */
2574 if (Sglext_isnotzero_low31p2(resultp2) ||
2575 Sglext_isone_lowp1(resultp1)) {
2576 /* either exactly half way and odd or
2577 * more than 1/2ulp */
2578 Sgl_increment(resultp1);
2579 }
2580 }
2581 break;
2582
2583 case ROUNDPLUS:
2584 if (Sgl_iszero_sign(resultp1)) {
2585 /* Round up positive results */
2586 Sgl_increment(resultp1);
2587 }
2588 break;
2589
2590 case ROUNDMINUS:
2591 if (Sgl_isone_sign(resultp1)) {
2592 /* Round down negative results */
2593 Sgl_increment(resultp1);
2594 }
2595
2596 case ROUNDZERO:;
2597 /* truncate is simple */
2598 } /* end switch... */
2599 if (Sgl_isone_hiddenoverflow(resultp1)) result_exponent++;
2600 }
2601 if (result_exponent >= SGL_INFINITY_EXPONENT) {
2602 /* Overflow */
2603 if (Is_overflowtrap_enabled()) {
2604 /*
2605 * Adjust bias of result
2606 */
2607 Sgl_setwrapped_exponent(resultp1,result_exponent,ovfl);
2608 Sgl_copytoptr(resultp1,dstptr);
2609 if (inexact)
2610 if (Is_inexacttrap_enabled())
2611 return (OPC_2E_OVERFLOWEXCEPTION |
2612 OPC_2E_INEXACTEXCEPTION);
2613 else Set_inexactflag();
2614 return (OPC_2E_OVERFLOWEXCEPTION);
2615 }
2616 inexact = TRUE;
2617 Set_overflowflag();
2618 Sgl_setoverflow(resultp1);
2619 } else if (result_exponent <= 0) { /* underflow case */
2620 if (Is_underflowtrap_enabled()) {
2621 /*
2622 * Adjust bias of result
2623 */
2624 Sgl_setwrapped_exponent(resultp1,result_exponent,unfl);
2625 Sgl_copytoptr(resultp1,dstptr);
2626 if (inexact)
2627 if (Is_inexacttrap_enabled())
2628 return (OPC_2E_UNDERFLOWEXCEPTION |
2629 OPC_2E_INEXACTEXCEPTION);
2630 else Set_inexactflag();
2631 return(OPC_2E_UNDERFLOWEXCEPTION);
2632 }
2633 else if (inexact && is_tiny) Set_underflowflag();
2634 }
2635 else Sgl_set_exponent(resultp1,result_exponent);
2636 Sgl_copytoptr(resultp1,dstptr);
2637 if (inexact)
2638 if (Is_inexacttrap_enabled()) return(OPC_2E_INEXACTEXCEPTION);
2639 else Set_inexactflag();
2640 return(NOEXCEPTION);
2641 }
2642
Linux with added FPC-III support