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