support.c
上传用户:nvosite88
上传日期:2007-01-17
资源大小:4983k
文件大小:9k
- /* 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);
- }
- }