support.c

<B>Digital的Unix操作系统VAX 4.2源码</B>

#ifndef lint
static char	*sccsid = "@(#)support.c	4.2	(ULTRIX)	12/4/90";
#endif lint

/************************************************************************
 *									*
 *			Copyright (c) 1986 by				*
 *		Digital Equipment Corporation, Maynard, MA		*
 *			All rights reserved.				*
 *									*
 *   This software is furnished under a license and may be used and	*
 *   copied  only  in accordance with the terms of such license and	*
 *   with the  inclusion  of  the  above  copyright  notice.   This	*
 *   software  or  any  other copies thereof may not be provided or	*
 *   otherwise made available to any other person.  No title to and	*
 *   ownership of the software is hereby transferred.			*
 *									*
 *   This software is  derived  from  software  received  from  the	*
 *   University    of   California,   Berkeley,   and   from   Bell	*
 *   Laboratories.  Use, duplication, or disclosure is  subject  to	*
 *   restrictions  under  license  agreements  with  University  of	*
 *   California and with AT&T.						*
 *									*
 *   The information in this software is subject to change  without	*
 *   notice  and should not be construed as a commitment by Digital	*
 *   Equipment Corporation.						*
 *									*
 *   Digital assumes no responsibility for the use  or  reliability	*
 *   of its software on equipment which is not supplied by Digital.	*
 *									*
 ************************************************************************/

/************************************************************************
*
*			Modification History
*
*		David Metsky,	17-Jan-86
*
* 001	Replaced old version with BSD 4.3 version as part of upgrade
*
*	Based on:	support.c 	1.1 (Berkeley)	5/23/85
*
*************************************************************************/

/* 
 * Some IEEE standard p754 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 p754 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 p754 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.
 */

#include <errno.h>
#include <math.h>

#ifdef VAX      /* VAX D format */
    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           /*IEEE double format */
    static unsigned short msign=0x7fff, mexp =0x7ff0  ;
    static short prep1=54, gap=4, bias=1025           ;
    static double novf=1.7E308, nunf=3.0E-308,zero=0.0;
#endif

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

#ifdef NATIONAL
        unsigned short *px=(unsigned short *) &x + 3;
#else /* VAX, SUN, ZILOG */
        unsigned short *px=(unsigned short *) &x;
#endif

        if( x == zero )  return(x); 

#ifdef VAX
        if( (k= *px & mexp ) != ~msign ) {
            if( N<-260) return(0.0); else if(N>260) return(HUGE_VAL);
#else   /* IEEE */
        if( (k= *px & mexp ) != mexp ) {
            if( N<-2100) return(0.0); else if(N>2100) return(HUGE_VAL);
            if( k == 0 ) {
                 x *= scalb(1.0,(int)prep1);  N -= prep1; return(scalb(x,N));}
#endif

            if((k = (k>>gap)+ N) > 0 )
                if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
                else {errno = EDOM; x=HUGE_VAL;}       /* overflow */
            else
                if( k > -prep1 ) 
                                        /* gradual underflow */
                    {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);}
                else
		    {errno = EDOM; return(0.0);}
            }
        return(x);
}


double copysign(x,y)
double x,y;
{
#ifdef NATIONAL
        unsigned short  *px=(unsigned short *) &x+3,
                        *py=(unsigned short *) &y+3;
#else /* VAX, SUN, ZILOG */
        unsigned short  *px=(unsigned short *) &x,
                        *py=(unsigned short *) &y;
#endif

/* Always guard against a reserved operand on vax. sin(-M_PI) cld IPO 04290 */
        if ( (*px & mexp) == 0 ) return(x);

        *px = ( *px & msign ) | ( *py & ~msign );
        return(x);
}

double logb(x)
double x; 
{

short *px=(short *) &x, k;

#ifdef VAX
        return( ((*px & mexp)>>gap) - bias);
#else /* IEEE */
        if( (k= *px & mexp ) != mexp )
            if ( k != 0 )
                return ( (k>>gap) - bias );
            else if( x != zero)
                return ( -1022.0 );
            else        
                return(-(0.0));    
#ifndef vax
        else if(x != x)
            return(x);
#endif
        else
            {*px &= msign; return(x);}
#endif
}

finite(x)
double x;    
{
return(1.0);
}

double drem(x,p)
double x,p;
{
        short sign;
        double hp,dp,tmp,drem(),scalb();
        unsigned short  k; 
#ifdef NATIONAL
        unsigned short
              *px=(unsigned short *) &x  +3, 
              *pp=(unsigned short *) &p  +3,
              *pd=(unsigned short *) &dp +3,
              *pt=(unsigned short *) &tmp+3;
#else /* VAX, SUN, ZILOG */
        unsigned short
              *px=(unsigned short *) &x  , 
              *pp=(unsigned short *) &p  ,
              *pd=(unsigned short *) &dp ,
              *pt=(unsigned short *) &tmp;
#endif

        *pp &= msign ;

#ifdef VAX
        if( ( *px & mexp ) == ~msign || p == zero )
#else /* IEEE */
        if( ( *px & mexp ) == mexp || p == zero )
#endif

#ifndef vax
                return( (x != x)? x:0.0 );
#else
		return(0.0);
#endif

        else  if ( ((*pp & mexp)>>gap) <= 1 ) 
                /* subnormal p, or almost subnormal p */
            { double b; b=scalb(1.0,(int)prep1);
              p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);}
        else  if ( p >= novf/2)
            { p /= 2 ; x /= 2; return(drem(x,p)*2);}
        else 
            {
                dp=p+p; hp=p/2;
                sign= *px & ~msign ;
                *px &= msign       ;
                while ( x > dp )
                    {
                        k=(*px & mexp) - (*pd & mexp) ;
                        tmp = dp ;
                        *pt += k ;

#ifdef VAX
                        if( x < tmp ) *pt -= 128 ;
#else /* IEEE */
                        if( x < tmp ) *pt -= 16 ;
#endif

                        x -= tmp ;
                    }
                if ( x > hp )
                    { x -= p ;  if ( x >= hp ) x -= p ; }

		*px = *px ^ sign;
                return( x);

            }
}