dfadd.c

这个linux源代码是很全面的~基本完整了~使用c编译的~由于时间问题我没有亲自测试~但就算用来做参考资料也是非常好的

/*
 * Linux/PA-RISC Project (http://www.parisc-linux.org/)
 *
 * Floating-point emulation code
 *  Copyright (C) 2001 Hewlett-Packard (Paul Bame) <[e-mail]>
 *
 *    This program is free software; you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation; either version 2, or (at your option)
 *    any later version.
 *
 *    This program is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with this program; if not, write to the Free Software
 *    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
/*
 * BEGIN_DESC
 *
 *  File:
 *	@(#)	pa/spmath/dfadd.c		$Revision: 1.1 $
 *
 *  Purpose:
 *	Double_add: add two double precision values.
 *
 *  External Interfaces:
 *	dbl_fadd(leftptr, rightptr, dstptr, status)
 *
 *  Internal Interfaces:
 *
 *  Theory:
 *	<<please update with a overview of the operation of this file>>
 *
 * END_DESC
*/


#include "float.h"
#include "dbl_float.h"

/*
 * Double_add: add two double precision values.
 */
dbl_fadd(
    dbl_floating_point *leftptr,
    dbl_floating_point *rightptr,
    dbl_floating_point *dstptr,
    unsigned int *status)
{
    register unsigned int signless_upper_left, signless_upper_right, save;
    register unsigned int leftp1, leftp2, rightp1, rightp2, extent;
    register unsigned int resultp1 = 0, resultp2 = 0;
    
    register int result_exponent, right_exponent, diff_exponent;
    register int sign_save, jumpsize;
    register boolean inexact = FALSE;
    register boolean underflowtrap;
        
    /* Create local copies of the numbers */
    Dbl_copyfromptr(leftptr,leftp1,leftp2);
    Dbl_copyfromptr(rightptr,rightp1,rightp2);

    /* A zero "save" helps discover equal operands (for later),  *
     * and is used in swapping operands (if needed).             */
    Dbl_xortointp1(leftp1,rightp1,/*to*/save);

    /*
     * check first operand for NaN's or infinity
     */
    if ((result_exponent = Dbl_exponent(leftp1)) == DBL_INFINITY_EXPONENT)
	{
	if (Dbl_iszero_mantissa(leftp1,leftp2)) 
	    {
	    if (Dbl_isnotnan(rightp1,rightp2)) 
		{
		if (Dbl_isinfinity(rightp1,rightp2) && save!=0) 
		    {
		    /* 
		     * invalid since operands are opposite signed infinity's
		     */
		    if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
                    Set_invalidflag();
                    Dbl_makequietnan(resultp1,resultp2);
		    Dbl_copytoptr(resultp1,resultp2,dstptr);
		    return(NOEXCEPTION);
		    }
		/*
	 	 * return infinity
	 	 */
		Dbl_copytoptr(leftp1,leftp2,dstptr);
		return(NOEXCEPTION);
		}
	    }
	else 
	    {
            /*
             * is NaN; signaling or quiet?
             */
            if (Dbl_isone_signaling(leftp1)) 
		{
               	/* trap if INVALIDTRAP enabled */
		if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
        	/* make NaN quiet */
        	Set_invalidflag();
        	Dbl_set_quiet(leftp1);
        	}
	    /* 
	     * is second operand a signaling NaN? 
	     */
	    else if (Dbl_is_signalingnan(rightp1)) 
		{
        	/* trap if INVALIDTRAP enabled */
               	if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
		/* make NaN quiet */
		Set_invalidflag();
		Dbl_set_quiet(rightp1);
		Dbl_copytoptr(rightp1,rightp2,dstptr);
		return(NOEXCEPTION);
		}
	    /*
 	     * return quiet NaN
 	     */
	    Dbl_copytoptr(leftp1,leftp2,dstptr);
 	    return(NOEXCEPTION);
	    }
	} /* End left NaN or Infinity processing */
    /*
     * check second operand for NaN's or infinity
     */
    if (Dbl_isinfinity_exponent(rightp1)) 
	{
	if (Dbl_iszero_mantissa(rightp1,rightp2)) 
	    {
	    /* return infinity */
	    Dbl_copytoptr(rightp1,rightp2,dstptr);
	    return(NOEXCEPTION);
	    }
        /*
         * is NaN; signaling or quiet?
         */
        if (Dbl_isone_signaling(rightp1)) 
	    {
            /* trap if INVALIDTRAP enabled */
	    if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
	    /* make NaN quiet */
	    Set_invalidflag();
	    Dbl_set_quiet(rightp1);
	    }
	/*
	 * return quiet NaN
 	 */
	Dbl_copytoptr(rightp1,rightp2,dstptr);
	return(NOEXCEPTION);
    	} /* End right NaN or Infinity processing */

    /* Invariant: Must be dealing with finite numbers */

    /* Compare operands by removing the sign */
    Dbl_copytoint_exponentmantissap1(leftp1,signless_upper_left);
    Dbl_copytoint_exponentmantissap1(rightp1,signless_upper_right);

    /* sign difference selects add or sub operation. */
    if(Dbl_ismagnitudeless(leftp2,rightp2,signless_upper_left,signless_upper_right))
	{
	/* Set the left operand to the larger one by XOR swap *
	 *  First finish the first word using "save"          */
	Dbl_xorfromintp1(save,rightp1,/*to*/rightp1);
	Dbl_xorfromintp1(save,leftp1,/*to*/leftp1);
     	Dbl_swap_lower(leftp2,rightp2);
	result_exponent = Dbl_exponent(leftp1);
	}
    /* Invariant:  left is not smaller than right. */ 

    if((right_exponent = Dbl_exponent(rightp1)) == 0)
        {
	/* Denormalized operands.  First look for zeroes */
	if(Dbl_iszero_mantissa(rightp1,rightp2)) 
	    {
	    /* right is zero */
	    if(Dbl_iszero_exponentmantissa(leftp1,leftp2))
		{
		/* Both operands are zeros */
		if(Is_rounding_mode(ROUNDMINUS))
		    {
		    Dbl_or_signs(leftp1,/*with*/rightp1);
		    }
		else
		    {
		    Dbl_and_signs(leftp1,/*with*/rightp1);
		    }
		}
	    else 
		{
		/* Left is not a zero and must be the result.  Trapped
		 * underflows are signaled if left is denormalized.  Result
		 * is always exact. */
		if( (result_exponent == 0) && Is_underflowtrap_enabled() )
		    {
		    /* need to normalize results mantissa */
	    	    sign_save = Dbl_signextendedsign(leftp1);
		    Dbl_leftshiftby1(leftp1,leftp2);
		    Dbl_normalize(leftp1,leftp2,result_exponent);
		    Dbl_set_sign(leftp1,/*using*/sign_save);
                    Dbl_setwrapped_exponent(leftp1,result_exponent,unfl);
		    Dbl_copytoptr(leftp1,leftp2,dstptr);
		    /* inexact = FALSE */
		    return(UNDERFLOWEXCEPTION);
		    }
		}
	    Dbl_copytoptr(leftp1,leftp2,dstptr);
	    return(NOEXCEPTION);
	    }

	/* Neither are zeroes */
	Dbl_clear_sign(rightp1);	/* Exponent is already cleared */
	if(result_exponent == 0 )
	    {
	    /* Both operands are denormalized.  The result must be exact
	     * and is simply calculated.  A sum could become normalized and a
	     * difference could cancel to a true zero. */
	    if( (/*signed*/int) save < 0 )
		{
		Dbl_subtract(leftp1,leftp2,/*minus*/rightp1,rightp2,
		/*into*/resultp1,resultp2);
		if(Dbl_iszero_mantissa(resultp1,resultp2))
		    {
		    if(Is_rounding_mode(ROUNDMINUS))
			{
			Dbl_setone_sign(resultp1);
			}
		    else
			{
			Dbl_setzero_sign(resultp1);
			}
		    Dbl_copytoptr(resultp1,resultp2,dstptr);
		    return(NOEXCEPTION);
		    }
		}
	    else
		{
		Dbl_addition(leftp1,leftp2,rightp1,rightp2,
		/*into*/resultp1,resultp2);
		if(Dbl_isone_hidden(resultp1))
		    {
		    Dbl_copytoptr(resultp1,resultp2,dstptr);
		    return(NOEXCEPTION);
		    }
		}
	    if(Is_underflowtrap_enabled())
		{
		/* need to normalize result */
	    	sign_save = Dbl_signextendedsign(resultp1);
		Dbl_leftshiftby1(resultp1,resultp2);
		Dbl_normalize(resultp1,resultp2,result_exponent);
		Dbl_set_sign(resultp1,/*using*/sign_save);
                Dbl_setwrapped_exponent(resultp1,result_exponent,unfl);
	        Dbl_copytoptr(resultp1,resultp2,dstptr);
		/* inexact = FALSE */
	        return(UNDERFLOWEXCEPTION);