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/dfsub.c		$Revision: 1.1 $
13 *
14 *  Purpose:
15 *	Double_subtract: subtract two double precision values.
16 *
17 *  External Interfaces:
18 *	dbl_fsub(leftptr, rightptr, dstptr, status)
19 *
20 *  Internal Interfaces:
21 *
22 *  Theory:
23 *	<<please update with a overview of the operation of this file>>
24 *
25 * END_DESC
26*/
27
28
29#include "float.h"
30#include "dbl_float.h"
31
32/*
33 * Double_subtract: subtract two double precision values.
34 */
35int
36dbl_fsub(
37	    dbl_floating_point *leftptr,
38	    dbl_floating_point *rightptr,
39	    dbl_floating_point *dstptr,
40	    unsigned int *status)
41    {
42    register unsigned int signless_upper_left, signless_upper_right, save;
43    register unsigned int leftp1, leftp2, rightp1, rightp2, extent;
44    register unsigned int resultp1 = 0, resultp2 = 0;
45
46    register int result_exponent, right_exponent, diff_exponent;
47    register int sign_save, jumpsize;
48    register boolean inexact = FALSE, underflowtrap;
49
50    /* Create local copies of the numbers */
51    Dbl_copyfromptr(leftptr,leftp1,leftp2);
52    Dbl_copyfromptr(rightptr,rightp1,rightp2);
53
54    /* A zero "save" helps discover equal operands (for later),  *
55     * and is used in swapping operands (if needed).             */
56    Dbl_xortointp1(leftp1,rightp1,/*to*/save);
57
58    /*
59     * check first operand for NaN's or infinity
60     */
61    if ((result_exponent = Dbl_exponent(leftp1)) == DBL_INFINITY_EXPONENT)
62	{
63	if (Dbl_iszero_mantissa(leftp1,leftp2))
64	    {
65	    if (Dbl_isnotnan(rightp1,rightp2))
66		{
67		if (Dbl_isinfinity(rightp1,rightp2) && save==0)
68		    {
69		    /*
70		     * invalid since operands are same signed infinity's
71		     */
72		    if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
73                    Set_invalidflag();
74                    Dbl_makequietnan(resultp1,resultp2);
75		    Dbl_copytoptr(resultp1,resultp2,dstptr);
76		    return(NOEXCEPTION);
77		    }
78		/*
79	 	 * return infinity
80	 	 */
81		Dbl_copytoptr(leftp1,leftp2,dstptr);
82		return(NOEXCEPTION);
83		}
84	    }
85	else
86	    {
87            /*
88             * is NaN; signaling or quiet?
89             */
90            if (Dbl_isone_signaling(leftp1))
91		{
92               	/* trap if INVALIDTRAP enabled */
93		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
94        	/* make NaN quiet */
95        	Set_invalidflag();
96        	Dbl_set_quiet(leftp1);
97        	}
98	    /*
99	     * is second operand a signaling NaN?
100	     */
101	    else if (Dbl_is_signalingnan(rightp1))
102		{
103        	/* trap if INVALIDTRAP enabled */
104               	if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
105		/* make NaN quiet */
106		Set_invalidflag();
107		Dbl_set_quiet(rightp1);
108		Dbl_copytoptr(rightp1,rightp2,dstptr);
109		return(NOEXCEPTION);
110		}
111	    /*
112 	     * return quiet NaN
113 	     */
114	    Dbl_copytoptr(leftp1,leftp2,dstptr);
115 	    return(NOEXCEPTION);
116	    }
117	} /* End left NaN or Infinity processing */
118    /*
119     * check second operand for NaN's or infinity
120     */
121    if (Dbl_isinfinity_exponent(rightp1))
122	{
123	if (Dbl_iszero_mantissa(rightp1,rightp2))
124	    {
125	    /* return infinity */
126	    Dbl_invert_sign(rightp1);
127	    Dbl_copytoptr(rightp1,rightp2,dstptr);
128	    return(NOEXCEPTION);
129	    }
130        /*
131         * is NaN; signaling or quiet?
132         */
133        if (Dbl_isone_signaling(rightp1))
134	    {
135            /* trap if INVALIDTRAP enabled */
136	    if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
137	    /* make NaN quiet */
138	    Set_invalidflag();
139	    Dbl_set_quiet(rightp1);
140	    }
141	/*
142	 * return quiet NaN
143 	 */
144	Dbl_copytoptr(rightp1,rightp2,dstptr);
145	return(NOEXCEPTION);
146    	} /* End right NaN or Infinity processing */
147
148    /* Invariant: Must be dealing with finite numbers */
149
150    /* Compare operands by removing the sign */
151    Dbl_copytoint_exponentmantissap1(leftp1,signless_upper_left);
152    Dbl_copytoint_exponentmantissap1(rightp1,signless_upper_right);
153
154    /* sign difference selects add or sub operation. */
155    if(Dbl_ismagnitudeless(leftp2,rightp2,signless_upper_left,signless_upper_right))
156	{
157	/* Set the left operand to the larger one by XOR swap *
158	 *  First finish the first word using "save"          */
159	Dbl_xorfromintp1(save,rightp1,/*to*/rightp1);
160	Dbl_xorfromintp1(save,leftp1,/*to*/leftp1);
161     	Dbl_swap_lower(leftp2,rightp2);
162	result_exponent = Dbl_exponent(leftp1);
163	Dbl_invert_sign(leftp1);
164	}
165    /* Invariant:  left is not smaller than right. */
166
167    if((right_exponent = Dbl_exponent(rightp1)) == 0)
168        {
169	/* Denormalized operands.  First look for zeroes */
170	if(Dbl_iszero_mantissa(rightp1,rightp2))
171	    {
172	    /* right is zero */
173	    if(Dbl_iszero_exponentmantissa(leftp1,leftp2))
174		{
175		/* Both operands are zeros */
176		Dbl_invert_sign(rightp1);
177		if(Is_rounding_mode(ROUNDMINUS))
178		    {
179		    Dbl_or_signs(leftp1,/*with*/rightp1);
180		    }
181		else
182		    {
183		    Dbl_and_signs(leftp1,/*with*/rightp1);
184		    }
185		}
186	    else
187		{
188		/* Left is not a zero and must be the result.  Trapped
189		 * underflows are signaled if left is denormalized.  Result
190		 * is always exact. */
191		if( (result_exponent == 0) && Is_underflowtrap_enabled() )
192		    {
193		    /* need to normalize results mantissa */
194	    	    sign_save = Dbl_signextendedsign(leftp1);
195		    Dbl_leftshiftby1(leftp1,leftp2);
196		    Dbl_normalize(leftp1,leftp2,result_exponent);
197		    Dbl_set_sign(leftp1,/*using*/sign_save);
198                    Dbl_setwrapped_exponent(leftp1,result_exponent,unfl);
199		    Dbl_copytoptr(leftp1,leftp2,dstptr);
200		    /* inexact = FALSE */
201		    return(UNDERFLOWEXCEPTION);
202		    }
203		}
204	    Dbl_copytoptr(leftp1,leftp2,dstptr);
205	    return(NOEXCEPTION);
206	    }
207
208	/* Neither are zeroes */
209	Dbl_clear_sign(rightp1);	/* Exponent is already cleared */
210	if(result_exponent == 0 )
211	    {
212	    /* Both operands are denormalized.  The result must be exact
213	     * and is simply calculated.  A sum could become normalized and a
214	     * difference could cancel to a true zero. */
215	    if( (/*signed*/int) save >= 0 )
216		{
217		Dbl_subtract(leftp1,leftp2,/*minus*/rightp1,rightp2,
218		 /*into*/resultp1,resultp2);
219		if(Dbl_iszero_mantissa(resultp1,resultp2))
220		    {
221		    if(Is_rounding_mode(ROUNDMINUS))
222			{
223			Dbl_setone_sign(resultp1);
224			}
225		    else
226			{
227			Dbl_setzero_sign(resultp1);
228			}
229		    Dbl_copytoptr(resultp1,resultp2,dstptr);
230		    return(NOEXCEPTION);
231		    }
232		}
233	    else
234		{
235		Dbl_addition(leftp1,leftp2,rightp1,rightp2,
236		 /*into*/resultp1,resultp2);
237		if(Dbl_isone_hidden(resultp1))
238		    {
239		    Dbl_copytoptr(resultp1,resultp2,dstptr);
240		    return(NOEXCEPTION);
241		    }
242		}
243	    if(Is_underflowtrap_enabled())
244		{
245		/* need to normalize result */
246	    	sign_save = Dbl_signextendedsign(resultp1);
247		Dbl_leftshiftby1(resultp1,resultp2);
248		Dbl_normalize(resultp1,resultp2,result_exponent);
249		Dbl_set_sign(resultp1,/*using*/sign_save);
250                Dbl_setwrapped_exponent(resultp1,result_exponent,unfl);
251		Dbl_copytoptr(resultp1,resultp2,dstptr);
252		/* inexact = FALSE */
253		return(UNDERFLOWEXCEPTION);
254		}
255	    Dbl_copytoptr(resultp1,resultp2,dstptr);
256	    return(NOEXCEPTION);
257	    }
258	right_exponent = 1;	/* Set exponent to reflect different bias
259				 * with denormalized numbers. */
260	}
261    else
262	{
263	Dbl_clear_signexponent_set_hidden(rightp1);
264	}
265    Dbl_clear_exponent_set_hidden(leftp1);
266    diff_exponent = result_exponent - right_exponent;
267
268    /*
269     * Special case alignment of operands that would force alignment
270     * beyond the extent of the extension.  A further optimization
271     * could special case this but only reduces the path length for this
272     * infrequent case.
273     */
274    if(diff_exponent > DBL_THRESHOLD)
275	{
276	diff_exponent = DBL_THRESHOLD;
277	}
278
279    /* Align right operand by shifting to right */
280    Dbl_right_align(/*operand*/rightp1,rightp2,/*shifted by*/diff_exponent,
281     /*and lower to*/extent);
282
283    /* Treat sum and difference of the operands separately. */
284    if( (/*signed*/int) save >= 0 )
285	{
286	/*
287	 * Difference of the two operands.  Their can be no overflow.  A
288	 * borrow can occur out of the hidden bit and force a post
289	 * normalization phase.
290	 */
291	Dbl_subtract_withextension(leftp1,leftp2,/*minus*/rightp1,rightp2,
292	 /*with*/extent,/*into*/resultp1,resultp2);
293	if(Dbl_iszero_hidden(resultp1))
294	    {
295	    /* Handle normalization */
296	    /* A straight forward algorithm would now shift the result
297	     * and extension left until the hidden bit becomes one.  Not
298	     * all of the extension bits need participate in the shift.
299	     * Only the two most significant bits (round and guard) are
300	     * needed.  If only a single shift is needed then the guard
301	     * bit becomes a significant low order bit and the extension
302	     * must participate in the rounding.  If more than a single
303	     * shift is needed, then all bits to the right of the guard
304	     * bit are zeros, and the guard bit may or may not be zero. */
305	    sign_save = Dbl_signextendedsign(resultp1);
306            Dbl_leftshiftby1_withextent(resultp1,resultp2,extent,resultp1,resultp2);
307
308            /* Need to check for a zero result.  The sign and exponent
309	     * fields have already been zeroed.  The more efficient test
310	     * of the full object can be used.
311	     */
312    	    if(Dbl_iszero(resultp1,resultp2))
313		/* Must have been "x-x" or "x+(-x)". */
314		{
315		if(Is_rounding_mode(ROUNDMINUS)) Dbl_setone_sign(resultp1);
316		Dbl_copytoptr(resultp1,resultp2,dstptr);
317		return(NOEXCEPTION);
318		}
319	    result_exponent--;
320	    /* Look to see if normalization is finished. */
321	    if(Dbl_isone_hidden(resultp1))
322		{
323		if(result_exponent==0)
324		    {
325		    /* Denormalized, exponent should be zero.  Left operand *
326		     * was normalized, so extent (guard, round) was zero    */
327		    goto underflow;
328		    }
329		else
330		    {
331		    /* No further normalization is needed. */
332		    Dbl_set_sign(resultp1,/*using*/sign_save);
333	    	    Ext_leftshiftby1(extent);
334		    goto round;
335		    }
336		}
337
338	    /* Check for denormalized, exponent should be zero.  Left    *
339	     * operand was normalized, so extent (guard, round) was zero */
340	    if(!(underflowtrap = Is_underflowtrap_enabled()) &&
341	       result_exponent==0) goto underflow;
342
343	    /* Shift extension to complete one bit of normalization and
344	     * update exponent. */
345	    Ext_leftshiftby1(extent);
346
347	    /* Discover first one bit to determine shift amount.  Use a
348	     * modified binary search.  We have already shifted the result
349	     * one position right and still not found a one so the remainder
350	     * of the extension must be zero and simplifies rounding. */
351	    /* Scan bytes */
352	    while(Dbl_iszero_hiddenhigh7mantissa(resultp1))
353		{
354		Dbl_leftshiftby8(resultp1,resultp2);
355		if((result_exponent -= 8) <= 0  && !underflowtrap)
356		    goto underflow;
357		}
358	    /* Now narrow it down to the nibble */
359	    if(Dbl_iszero_hiddenhigh3mantissa(resultp1))
360		{
361		/* The lower nibble contains the normalizing one */
362		Dbl_leftshiftby4(resultp1,resultp2);
363		if((result_exponent -= 4) <= 0 && !underflowtrap)
364		    goto underflow;
365		}
366	    /* Select case were first bit is set (already normalized)
367	     * otherwise select the proper shift. */
368	    if((jumpsize = Dbl_hiddenhigh3mantissa(resultp1)) > 7)
369		{
370		/* Already normalized */
371		if(result_exponent <= 0) goto underflow;
372		Dbl_set_sign(resultp1,/*using*/sign_save);
373		Dbl_set_exponent(resultp1,/*using*/result_exponent);
374		Dbl_copytoptr(resultp1,resultp2,dstptr);
375		return(NOEXCEPTION);
376		}
377	    Dbl_sethigh4bits(resultp1,/*using*/sign_save);
378	    switch(jumpsize)
379		{
380		case 1:
381		    {
382		    Dbl_leftshiftby3(resultp1,resultp2);
383		    result_exponent -= 3;
384		    break;
385		    }
386		case 2:
387		case 3:
388		    {
389		    Dbl_leftshiftby2(resultp1,resultp2);
390		    result_exponent -= 2;
391		    break;
392		    }
393		case 4:
394		case 5:
395		case 6:
396		case 7:
397		    {
398		    Dbl_leftshiftby1(resultp1,resultp2);
399		    result_exponent -= 1;
400		    break;
401		    }
402		}
403	    if(result_exponent > 0)
404		{
405		Dbl_set_exponent(resultp1,/*using*/result_exponent);
406		Dbl_copytoptr(resultp1,resultp2,dstptr);
407		return(NOEXCEPTION);		/* Sign bit is already set */
408		}
409	    /* Fixup potential underflows */
410	  underflow:
411	    if(Is_underflowtrap_enabled())
412		{
413		Dbl_set_sign(resultp1,sign_save);
414                Dbl_setwrapped_exponent(resultp1,result_exponent,unfl);
415		Dbl_copytoptr(resultp1,resultp2,dstptr);
416		/* inexact = FALSE */
417		return(UNDERFLOWEXCEPTION);
418		}
419	    /*
420	     * Since we cannot get an inexact denormalized result,
421	     * we can now return.
422	     */
423	    Dbl_fix_overshift(resultp1,resultp2,(1-result_exponent),extent);
424	    Dbl_clear_signexponent(resultp1);
425	    Dbl_set_sign(resultp1,sign_save);
426	    Dbl_copytoptr(resultp1,resultp2,dstptr);
427	    return(NOEXCEPTION);
428	    } /* end if(hidden...)... */
429	/* Fall through and round */
430	} /* end if(save >= 0)... */
431    else
432	{
433	/* Subtract magnitudes */
434	Dbl_addition(leftp1,leftp2,rightp1,rightp2,/*to*/resultp1,resultp2);
435	if(Dbl_isone_hiddenoverflow(resultp1))
436	    {
437	    /* Prenormalization required. */
438	    Dbl_rightshiftby1_withextent(resultp2,extent,extent);
439	    Dbl_arithrightshiftby1(resultp1,resultp2);
440	    result_exponent++;
441	    } /* end if hiddenoverflow... */
442	} /* end else ...subtract magnitudes... */
443
444    /* Round the result.  If the extension is all zeros,then the result is
445     * exact.  Otherwise round in the correct direction.  No underflow is
446     * possible. If a postnormalization is necessary, then the mantissa is
447     * all zeros so no shift is needed. */
448  round:
449    if(Ext_isnotzero(extent))
450	{
451	inexact = TRUE;
452	switch(Rounding_mode())
453	    {
454	    case ROUNDNEAREST: /* The default. */
455	    if(Ext_isone_sign(extent))
456		{
457		/* at least 1/2 ulp */
458		if(Ext_isnotzero_lower(extent)  ||
459		  Dbl_isone_lowmantissap2(resultp2))
460		    {
461		    /* either exactly half way and odd or more than 1/2ulp */
462		    Dbl_increment(resultp1,resultp2);
463		    }
464		}
465	    break;
466
467	    case ROUNDPLUS:
468	    if(Dbl_iszero_sign(resultp1))
469		{
470		/* Round up positive results */
471		Dbl_increment(resultp1,resultp2);
472		}
473	    break;
474
475	    case ROUNDMINUS:
476	    if(Dbl_isone_sign(resultp1))
477		{
478		/* Round down negative results */
479		Dbl_increment(resultp1,resultp2);
480		}
481
482	    case ROUNDZERO:;
483	    /* truncate is simple */
484	    } /* end switch... */
485	if(Dbl_isone_hiddenoverflow(resultp1)) result_exponent++;
486	}
487    if(result_exponent == DBL_INFINITY_EXPONENT)
488        {
489        /* Overflow */
490        if(Is_overflowtrap_enabled())
491	    {
492	    Dbl_setwrapped_exponent(resultp1,result_exponent,ovfl);
493	    Dbl_copytoptr(resultp1,resultp2,dstptr);
494	    if (inexact)
495	    if (Is_inexacttrap_enabled())
496		return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
497		else Set_inexactflag();
498	    return(OVERFLOWEXCEPTION);
499	    }
500        else
501	    {
502	    inexact = TRUE;
503	    Set_overflowflag();
504	    Dbl_setoverflow(resultp1,resultp2);
505	    }
506	}
507    else Dbl_set_exponent(resultp1,result_exponent);
508    Dbl_copytoptr(resultp1,resultp2,dstptr);
509    if(inexact)
510	if(Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
511	else Set_inexactflag();
512    return(NOEXCEPTION);
513    }
514