ansimath.c

VXWORKS源代码

static long        p2x[] = { _0x(0b60,bc36), _0x(ec94,b5f5)};
static long        p3x[] = { _0x(b355,398a), _0x(f15f,792e)};
static long        p4x[] = { _0x(ea0e,b6dd), _0x(5f84,2e93)};
static long        p5x[] = { _0x(bb4b,3431), _0x(2683,95f5)};
#define    ln2hi    (*(double*)ln2hix)
#define    ln2lo    (*(double*)ln2lox)
#define   lnhuge    (*(double*)lnhugex)
#define   lntiny    (*(double*)lntinyx)
#define   invln2    (*(double*)invln2x)
#define       p1    (*(double*)p1x)
#define       p2    (*(double*)p2x)
#define       p3    (*(double*)p3x)
#define       p4    (*(double*)p4x)
#define       p5    (*(double*)p5x)

#else	/* defined(vax)||defined(tahoe) */
static double
p1     =  1.6666666666666601904E-1    , /*Hex  2^-3    *  1.555555555553E */
p2     = -2.7777777777015593384E-3    , /*Hex  2^-9    * -1.6C16C16BEBD93 */
p3     =  6.6137563214379343612E-5    , /*Hex  2^-14   *  1.1566AAF25DE2C */
p4     = -1.6533902205465251539E-6    , /*Hex  2^-20   * -1.BBD41C5D26BF1 */
p5     =  4.1381367970572384604E-8    , /*Hex  2^-25   *  1.6376972BEA4D0 */
ln2hi  =  6.9314718036912381649E-1    , /*Hex  2^ -1   *  1.62E42FEE00000 */
ln2lo  =  1.9082149292705877000E-10   , /*Hex  2^-33   *  1.A39EF35793C76 */
lnhuge =  7.1602103751842355450E2     , /*Hex  2^  9   *  1.6602B15B7ECF2 */
lntiny = -7.5137154372698068983E2     , /*Hex  2^  9   * -1.77AF8EBEAE354 */
invln2 =  1.4426950408889633870E0     ; /*Hex  2^  0   *  1.71547652B82FE */
#endif	/* defined(vax)||defined(tahoe) */

/*****************************************************************************
*
* exp - compute an exponential value (ANSI)
*
* This routine returns the exponential value of <x> in
* double precision (IEEE double, 53 bits).
*
* A range error occurs if <x> is too large.
*
* INTERNAL:
* Method:
* (1) Argument Reduction: given the input <x>, find <r> and integer <k>
*     such that:
*
*         x = k*ln2 + r,  |r| <= 0.5*ln2
* 
*     <r> will be represented as r := z+c for better accuracy.
* 
* (2) Compute exp(r) by
*
*         exp(r) = 1 + r + r*R1/(2-R1)
*
*     where:
*
*         R1 = x - x^2*(p1+x^2*(p2+x^2*(p3+x^2*(p4+p5*x^2))))
*
* (3)     exp(x) = 2^k * exp(r)
*
* INCLUDE FILES: math.h
*
* RETURNS: The double-precision exponential value of <x>.
*
* Special cases:
*     If <x> is +INF or NaN, exp() returns <x>.
*     If <x> is -INF, it returns 0.
*
* SEE ALSO: mathALib
*
* INTERNAL:
* Coded in C by K.C. Ng, 1/19/85;
* Revised by K.C. Ng on 2/6/85, 2/15/85, 3/7/85, 3/24/85, 4/16/85, 6/14/86.
*/

double exp
    (
    double x	/* exponent */
    )

    {
	double scalb(), copysign(), z,hi,lo,c;
	int k,finite();

#if !defined(vax)&&!defined(tahoe)
	if(x!=x) return(x);	/* x is NaN */
#endif	/* !defined(vax)&&!defined(tahoe) */
	if( x <= lnhuge ) {
		if( x >= lntiny ) {

		    /* argument reduction : x --> x - k*ln2 */

			k=invln2*x+copysign(0.5,x);	/* k=NINT(x/ln2) */

		    /* express x-k*ln2 as hi-lo and let x=hi-lo rounded */

			hi=x-k*ln2hi;
			x=hi-(lo=k*ln2lo);

		    /* return 2^k*[1+x+x*c/(2+c)]  */
			z=x*x;
			c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5))));
			return  scalb(1.0+(hi-(lo-(x*c)/(2.0-c))),k);

		}
		/* end of x > lntiny */

		else
		     /* exp(-big#) underflows to zero */
		     if(finite(x))  return(scalb(1.0,-5000));

		     /* exp(-INF) is zero */
		     else return(0.0);
	}
	/* end of x < lnhuge */

	else
	/* exp(INF) is INF, exp(+big#) overflows to INF */
	    return( finite(x) ?  scalb(1.0,5000)  : x);
    }
/* fabs.c - fabs math routine */

/* Copyright 1992-1993 Wind River Systems, Inc. */

/*
modification history
--------------------
01f,05feb93,jdi  doc changes based on kdl review.
01e,02dec92,jdi  doc tweaks.
01d,28oct92,jdi  documentation cleanup.
01c,20sep92,smb  documentation additions
01b,30jul92,kdl  changed _d_type() calls to fpTypeGet().
01a,08jul92,smb  documentation.
*/

/*
DESCRIPTION
SEE ALSO: American National Standard X3.159-1989
NOMANUAL
*/

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


/*******************************************************************************
*
* fabs - compute an absolute value (ANSI)
*
* This routine returns the absolute value of <v> in double precision.
*
* INCLUDE FILES: math.h
*
* RETURNS: The double-precision absolute value of <v>.
*
* ERRNO: EDOM, ERANGE
*
* SEE ALSO: mathALib
*/

double fabs
    (
    double v	/* number to return the absolute value of */
    )

    {
    switch (fpTypeGet (v, NULL))
	{
	case ZERO:			/* ZERO */
	    return (0.0);

	case NAN:			/* NOT  A NUMBER */
	    errno = EDOM;
	    return (v);

	case INF:			/* INFINITE */
	    errno = ERANGE;

	case INTEGER:			/* INTEGER */
	case REAL:			/* REAL */
	    return ((v < 0.0) ? - v : v);

	default:
	    return (0.0);
	}
    }
/* floor.c - floor math routine */

/* Copyright 1992-1993 Wind River Systems, Inc. */

/*
modification history
--------------------
01f,05feb93,jdi  doc changes based on kdl review.
01e,02dec92,jdi  doc tweaks.
01d,28oct92,jdi  documentation cleanup.
01c,20sep92,smb  documentation additions
01b,30jul92,kdl  changed _d_type() calls to fpTypeGet().
01a,08jul92,smb  documentation.
*/

/*
DESCRIPTION
SEE ALSO: American National Standard X3.159-1989
NOMANUAL
*/

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


/*******************************************************************************
*
* floor - compute the largest integer less than or equal to a specified value (ANSI)
*
* This routine returns the largest integer less than or equal to <v>, in
* double precision.
* 
* INCLUDE FILES: math.h
*
* RETURNS:
* The largest integral value less than or equal to <v>, in double precision.
*
* SEE ALSO: mathALib
*/

double floor
    (
    double v	/* value to find the floor of */
    )

    {
    double r;

    switch (fpTypeGet (v, &r))	/* find out the type of v */
	{
	case ZERO:
	case INF:	return (0);

	case NAN:	return (v);

	case INTEGER:	return (v);

	case REAL:	return (((v < 0.0) && (v != r)) ? (r - 1.0) : (r));

	default:	return (0);		/* this should never happen */
	}
    }
/* fmod.c - fmod math routine */

/* Copyright 1992-1993 Wind River Systems, Inc. */

/*
modification history
--------------------
01f,05feb93,jdi  doc changes based on kdl review.
01e,02dec92,jdi  doc tweaks.
01d,28oct92,jdi  documentation cleanup.
01c,20sep92,smb  documentation additions
01b,30jul92,kdl  changed _d_type() calls to fpTypeGet().
01a,08jul92,smb  documentation.
*/

/*
DESCRIPTION
SEE ALSO: American National Standard X3.159-1989
NOMANUAL
*/

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


/*******************************************************************************
*
* fmod - compute the remainder of x/y (ANSI)
*
* This routine returns the remainder of <x>/<y> with the sign of <x>,
* in double precision.
* 
* INCLUDE FILES: math.h
*
* RETURNS: The value <x> - <i> * <y>, for some integer <i>.  If <y> is
* non-zero, the result has the same sign as <x> and magnitude less than the
* magnitude of <y>.  If <y> is zero, fmod() returns zero.
*
* ERRNO: EDOM
*
* SEE ALSO: mathALib
*/

double fmod
    (
    double x,	/* numerator   */
    double y	/* denominator */
    )

    {
    double	t;
    short	negative = 0;
    int		errx = fpTypeGet (x, NULL);	/* determine number type */
    int		erry = fpTypeGet (y, NULL);	/* determine number type */

    if (errx == NAN || erry == NAN || errx == INF || erry == ZERO)
        {
        errno = EDOM;				/* Domain error */
        return ((errx == NAN) ? (x) : ((erry == NAN) ? (y) : (0)));
        }
    
    if (errx == ZERO || erry == INF)
	return (x);		

    /* make x and y absolute */

    if (y < 0.0)
	y = -y;

    if (x < 0.0)
	{
	x = -x;
	negative = 1;
	}

    /* loop substracting y from x until a value less than y remains */

    for (t = x; t > y; t -= y)
        ;

    return ((t == y) ? (0.0) : ((negative) ? (-t) : (t)));
    }
/* frexp.c - frexp math routine */

/* Copyright 1992-1996 Wind River Systems, Inc. */

/*
modification history
--------------------
01h,05aug96,dbt  frexp now returned a value >= 0.5 and strictly less than 1.0.
		 (SPR #4280).
		 Updated copyright.
01g,05feb93,jdi  doc changes based on kdl review.
01f,02dec92,jdi  doc tweaks.
01e,28oct92,jdi  documentation cleanup.
01d,21sep92,smb  replaced orignal versions.
01c,20sep92,smb  documentation additions
01b,30jul92,kdl  replaced routine contents with frexp() from floatLib.c.
01a,08jul92,smb  written, documentation.
*/

/*
DESCRIPTION

SEE ALSO: American National Standard X3.159-1989

NOMANUAL
*/

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

/*******************************************************************************
*
* frexp - break a floating-point number into a normalized fraction and power of 2 (ANSI)
*
* This routine breaks a double-precision number <value> into a normalized
* fraction and integral power of 2.  It stores the integer exponent in <pexp>.
*
* INCLUDE FILES: math.h
*
* RETURNS: 
* The double-precision value <x>, such that the magnitude of <x> is
* in the interval [1/2,1) or zero, and <value> equals <x> times
* 2 to the power of <pexp>. If <value> is zero, both parts of the result
* are zero.
*
* ERRNO: EDOM
*
*/
double frexp
    (
    double value,	/* number to be normalized */
    int *pexp   	/* pointer to the exponent */
    )

    {
    double r;

    *pexp = 0;

    /* determine number type */

    switch (fpTypeGet(value, &r))	
      	{
        case NAN:		/* NOT A NUMBER */
    	case INF:		/* INFINITE */
            errno = EDOM;
	    return (value);

    	case ZERO:		/* ZERO */
	    return (0.0);

	default:
	    /*
	     * return value must be strictly less than 1.0 and >=0.5 .
	     */
	    if ((r = fabs(value)) >= 1.0)
	        for(; (r >= 1.0); (*pexp)++, r /= 2.0);
	    else
	        for(; (r < 0.5); (*pexp)--, r *= 2.0);

	    return (value < 0 ? -r : r);
        }
    }
/* ldexp.c - ldexp math routine */

/* Copyright 1992-1993 Wind River Systems, Inc. */

/*
modification history
--------------------
01g,05feb93,jdi  doc changes based on kdl review.
01f,02dec92,jdi  doc tweaks.
01e,28oct92,jdi  documentation cleanup.
01d,21sep92,smb  replaced original version.
01c,20sep92,smb  documentation additions
01b,30jul92,kdl  replaced routine contents with ldexp() from floatLib.c.
01a,08jul92,smb  written, documentation.
*/

/*
DESCRIPTION

SEE ALSO: American National Standard X3.159-1989