cachepc-linux

Fork of AMDESE/linux with modifications for CachePC side-channel attack
git clone https://git.sinitax.com/sinitax/cachepc-linux
Log | Files | Refs | README | LICENSE | sfeed.txt

dfadd.c (15114B)


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