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