support.c

VXWORKS源代码

/* support.c - math routines */

/* Copyright 1991-2001 Wind River Systems, Inc. */

/*
modification history
--------------------
01j,03oct01,to   use IEEE little endian for ARM
01i,20apr01,h_k  fixed P_INDEX for ARM big_endian.
01h,11jan99,dra  removed non-ANSI aliasing code.
01g,12feb97,jpd  added support for ARM with early Cygnus release.
01f,23sep92,kdl  removed define of sparcHardSqrt (now in sqrt.c).
01e,23sep92,kdl  removed sqrt(), placed in sqrt.c.
01d,20sep92,smb  documentation additions.
01c,08jul92,smb  merged with ANSI and documentation.
                 set EDOM in sqrt().
01b,02jul92,jwt  changed sysHardSqrt to sparcHardSqrt; removed warnings.
01a,25jun91,jwt  created sqrt() function select - hard versus soft.
*/

/*
DESCRIPTION
* Copyright (c) 1985 Regents of the University of California.
* All rights reserved.
*
* Redistribution and use in source and binary forms are permitted
* provided that the above copyright notice and this paragraph are
* duplicated in all such forms and that any documentation,
* advertising materials, and other materials related to such
* distribution and use acknowledge that the software was developed
* by the University of California, Berkeley.  The name of the
* University may not be used to endorse or promote products derived
* from this software without specific prior written permission.
* THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
* WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
*
* All recipients should regard themselves as participants in an ongoing
* research project and hence should feel obligated to report their
* experiences (good or bad) with these elementary function codes, using
* the sendbug(8) program, to the authors.
*
* Some IEEE standard 754 recommended functions and remainder and sqrt for
* supporting the C elementary functions.
******************************************************************************
* WARNING:
*      These codes are developed (in double) to support the C elementary
* functions temporarily. They are not universal, and some of them are very
* slow (in particular, drem and sqrt is extremely inefficient). Each
* computer system should have its implementation of these functions using
* its own assembler.
******************************************************************************
*
* IEEE 754 required operations:
*     drem(x,p)
*              returns  x REM y  =  x - [x/y]*y , where [x/y] is the integer
*              nearest x/y; in half way case, choose the even one.
*     sqrt(x)
*              returns the square root of x correctly rounded according to
*		the rounding mod.
*
* IEEE 754 recommended functions:
* (a) copysign(x,y)
*              returns x with the sign of y.
* (b) scalb(x,N)
*              returns  x * (2**N), for integer values N.
* (c) logb(x)
*              returns the unbiased exponent of x, a signed integer in
*              double precision, except that logb(0) is -INF, logb(INF)
*              is +INF, and logb(NAN) is that NAN.
* (d) finite(x)
*              returns the value TRUE if -INF < x < +INF and returns
*              FALSE otherwise.
*
*
* CODED IN C BY K.C. NG, 11/25/84;
* REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
*
* 25jun91,jwt  created sqrt() function select - hard versus soft
SEE ALSO: American National Standard X3.159-1989
NOMANUAL
*/

#include "vxWorks.h"
#include "math.h"
#include "private/mathP.h"
#include "errno.h"

#if (_BYTE_ORDER == _LITTLE_ENDIAN)
#define P_INDEX	3
#else	/* little_endian */
#define P_INDEX	0
#endif	/* little_endian */

typedef union
    {
    unsigned short	s[4];
    double		d;
    } tdouble;

#if defined(vax)||defined(tahoe)      /* VAX D format */
#include <errno.h>
    static unsigned short msign=0x7fff , mexp =0x7f80 ;
    static short  prep1=57, gap=7, bias=129           ;
    static double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
#else	/* defined(vax)||defined(tahoe) */
    static unsigned short msign=0x7fff, mexp =0x7ff0  ;
    static short prep1=54, gap=4, bias=1023           ;
    static double novf=1.7E308, nunf=3.0E-308,zero=0.0;
#endif	/* defined(vax)||defined(tahoe) */

/*****************************************************************************
* scalb	-
*
* RETURN:
* NOMANUAL
*/

double scalb(x,N)
double x; int N;
{
        int k;
        double scalb();
	tdouble y;

	y.d = x;
        if( y.d == zero )
	    return(y.d);

#if defined(vax)||defined(tahoe)
        if( (k= y.s[P_INDEX] & mexp ) != ~msign )
	    {
            if (N < -260)
		return(nunf*nunf);
	    else if (N > 260) {
		extern double infnan(),copysign();
		return(copysign(infnan(ERANGE),y.d));
	    }
#else	/* defined(vax)||defined(tahoe) */
        if( (k= y.s[P_INDEX] & mexp ) != mexp )
	    {
            if( N<-2100)
		return(nunf*nunf);
	    else if(N>2100)
		return(novf+novf);
            if( k == 0 )
		{
		y.d *= scalb(1.0,(int)prep1);
		N -= prep1;
		return(scalb(y.d,N));
		}
#endif	/* defined(vax)||defined(tahoe) */

            if((k = (k>>gap)+ N) > 0 )
                if( k < (mexp>>gap) )
		    y.s[P_INDEX] = (y.s[P_INDEX]&~mexp) | (k<<gap);
                else
		    x=novf+novf;               /* overflow */
            else
                if( k > -prep1 )
                                        /* gradual underflow */
                    {
		    y.s[P_INDEX]=(y.s[P_INDEX]&~mexp)|(short)(1<<gap);
		    y.d *= scalb(1.0,k-1);
		    }
                else
                return(nunf*nunf);
            }
        return(y.d);
}


/*****************************************************************************
* copysign - determine the sign
*
* RETURN:
* NOMANUAL
*/

double copysign(x,y)
double x,y;
{
	tdouble a, b;

	a.d = x;
	b.d = y;

#if defined(vax)||defined(tahoe)
        if ( (a.s[P_INDEX] & mexp) == 0 ) return(a.d);
#endif	/* defined(vax)||defined(tahoe) */

        a.s[P_INDEX] = ( a.s[P_INDEX] & msign ) | ( b.s[P_INDEX] & ~msign );
        return(a.d);
}

/*****************************************************************************
* logb -
*
* RETURN:
* NOMANUAL
*/

double logb(x)
double x;
{
	tdouble a;
	short k;

	a.d = x;
#if defined(vax)||defined(tahoe)
        return (int)(((a.s[P_INDEX]&mexp)>>gap)-bias);
#else	/* defined(vax)||defined(tahoe) */
        if( (k= a.s[P_INDEX] & mexp ) != mexp )
            if ( k != 0 )
                return ( (k>>gap) - bias );
            else if( a.d != zero)
                return ( -1022.0 );
            else
                return(-(1.0/zero));
        else if(a.d != a.d)
            return(a.d);
        else
            {a.s[P_INDEX] &= msign; return(a.d);}
#endif	/* defined(vax)||defined(tahoe) */
}

/*****************************************************************************
* finite - the finite value for this machine 
*
* RETURN:
* NOMANUAL
*/

int finite(x)
double x;
{
	tdouble a;
	a.d = x;
#if defined(vax)||defined(tahoe)
        return(1);
#else	/* defined(vax)||defined(tahoe) */
        return( ( a.s[P_INDEX] & mexp ) != mexp );
#endif	/* defined(vax)||defined(tahoe) */
}

/*****************************************************************************
* drem - remainder
*
* RETURN:
* NOMANUAL
*/

double drem(x,p)
double x,p;
{
        short sign;
        double hp, drem(),scalb();
	tdouble tmp, xd, pd, dp;
        unsigned short  k;

	xd.d = x;
	pd.d = p;

        pd.s[P_INDEX] &= msign;

#if defined(vax)||defined(tahoe)
        if( ( xd.s[P_INDEX] & mexp ) == ~msign )	/* is x a reserved operand? */
#else	/* defined(vax)||defined(tahoe) */
        if( ( xd.s[P_INDEX] & mexp ) == mexp )
#endif	/* defined(vax)||defined(tahoe) */
		return  (xd.d-pd.d)-(xd.d-pd.d);	/* create nan if x is inf */
	if (pd.d == zero) {
#if defined(vax)||defined(tahoe)
		extern double infnan();
		return(infnan(EDOM));
#else	/* defined(vax)||defined(tahoe) */
		return zero/zero;
#endif	/* defined(vax)||defined(tahoe) */
	}

#if defined(vax)||defined(tahoe)
        if( ( pd.s[P_INDEX] & mexp ) == ~msign )	/* is p a reserved operand? */
#else	/* defined(vax)||defined(tahoe) */
        if( ( pd.s[P_INDEX] & mexp ) == mexp )
#endif	/* defined(vax)||defined(tahoe) */
		{ if (pd.d != pd.d) return pd.d; else return xd.d;}

        else  if ( ((pd.s[P_INDEX] & mexp)>>gap) <= 1 )
                /* subnormal p, or almost subnormal p */
            { double b; b=scalb(1.0,(int)prep1);
              pd.d *= b; xd.d = drem(xd.d,pd.d); xd.d *= b; return(drem(xd.d,pd.d)/b);}
        else  if ( pd.d >= novf/2)
            { pd.d /= 2 ; xd.d /= 2; return(drem(xd.d,pd.d)*2);}
        else
            {
                dp.d=pd.d+pd.d;
		hp=pd.d/2;
                sign= xd.s[P_INDEX] & ~msign ;
                xd.s[P_INDEX] &= msign       ;
                while ( xd.d > dp.d )
                    {
                        k=(xd.s[P_INDEX] & mexp) - (dp.s[P_INDEX] & mexp) ;
                        tmp.d = dp.d ;
                        tmp.s[P_INDEX] += k ;

#if defined(vax)||defined(tahoe)
                        if( xd.d < tmp.d ) tmp.s[P_INDEX] -= 128 ;
#else	/* defined(vax)||defined(tahoe) */
                        if( xd.d < tmp.d ) tmp.s[P_INDEX] -= 16 ;
#endif	/* defined(vax)||defined(tahoe) */

                        xd.d -= tmp.d ;
                    }
                if ( xd.d > hp )
                    { xd.d -= pd.d ;  if ( xd.d >= hp ) xd.d -= pd.d ; }

#if defined(vax)||defined(tahoe)
		if (xd.d)
#endif	/* defined(vax)||defined(tahoe) */
			xd.s[P_INDEX] ^= sign;
                return( xd.d);

            }
}