ARM: 6979/1: mach-vt8500: add forgotten irq_data conversion
[linux-2.6.git] / arch / parisc / math-emu / fmpyfadd.c
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) <bame@debian.org>
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