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