mathsoftlib.c

vxworks source code, used for develop vxworks system.

    if (dblParam < 0.0)
	return (result + PI);		/* [pi/2 .. pi] */
    else
	return (result);		/* [0 .. pi/2) */
    }

/******************************************************************************
*
* mathSoftAtan2 - software floating point arc tangent of two arguments
*
* This routine takes two input double-precision floating point
* parameters, <dblY> and <dblX>, and returns the double precision
* arc tangent of <dblY> / <dblX>.
*
* RETURNS: double-precision arc tangent of <dblY> / <dblX>.
*/

LOCAL double mathSoftAtan2 
    (
    double	dblY,
    double	dblX 
    )
    {

    double	result;


    result = atan (dblY/dblX);

    if (dblX >= 0.0)
	return (result);		/* [-pi/2 .. pi/2) */
    else if (result < 0.0)
	return (result + PI);		/* [pi/2 .. pi) */
    else
	return (result - PI);		/* [-pi .. -pi/2) */
    }

/******************************************************************************
*
* mathSoftHypot - software floating point Euclidean distance
*
* This routine takes two input double-precision floating point
* parameters and returns length of the corresponding Euclidean distance
* (hypotenuse).
*
* RETURNS: double-precision hypotenuse.
*/

LOCAL double mathSoftHypot 
    (
    double	dblX,
    double	dblY 
    )
    {

    return (sqrt (dblX * dblX  +  dblY * dblY));
    }


/******************************************************************************
*
* mathSoftCbrt - software floating point cube root
*
* This routine takes the input double-precision floating point
* parameter and returns the cube root.
*
* kahan's cube root (53 bits IEEE double precision)
* for IEEE machines only
* coded in C by K.C. Ng, 4/30/85
*
* Accuracy:
*	better than 0.667 ulps according to an error analysis. Maximum
* error observed was 0.666 ulps in an 1,000,000 random arguments test.
*
* Warning: this code is semi machine dependent; the ordering of words in
* a floating point number must be known in advance. I assume that the
* long interger at the address of a floating point number will be the
* leading 32 bits of that floating point number (i.e., sign, exponent,
* and the 20 most significant bits).
* On a National machine, it has different ordering; therefore, this code
* must be compiled with flag -DNATIONAL.
*
* RETURNS: double-precision cube root.
*
* AUTHOR:  Kahan's cube root, coded in C by K.C. Ng, UC Berkeley, 4/30/85.
*	   Copyright (c) 1985 Regents of the University of California.
*	   Adapted by Scott Huddleston, Computer Research Lab, Tektronix.
*
* Use and reproduction of this software are granted  in  accordance  with
* the terms and conditions specified in  the  Berkeley  Software  License
* Agreement (in particular, this entails acknowledgement of the programs'
* source, and inclusion of this notice) with the additional understanding
* that  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 "sendbug 4bsd-bugs@BERKELEY", to the authors.
*/

LOCAL double mathSoftCbrt 
    (
    double	x 
    )
    {
	double r,s,t=0.0,w;
	unsigned long *px = (unsigned long *) &x,
	              *pt = (unsigned long *) &t,
		      mexp,sign;

static long B0 = 0x43500000;	/* (double) 2**54 */
static long B1 = 0x2a9f7893,
			B2 = 0x297f7893;	/* B1 / (double) 2**18 */

static double
	    C= 19./35.,
	    D= -864./1225.,
	    E= 99./70.,
	    F= 45./28.,
	    G= 5./14.;

#ifdef NATIONAL /* ordering of words in a floating points number */
	int n0=1,n1=0;
#else
	int n0=0,n1=1;
#endif

	mexp = px[n0]&0x7ff00000;
	if(isINFNaN(x)) return(x); /* cbrt(NaN,INF) is itself */
	if(x==0.0) return(x);	/* cbrt(0) is itself */

	sign=px[n0]&0x80000000; /* sign= sign(x) */
	px[n0] ^= sign;		/* x=|x| */

    /* rough cbrt to 5 bits */
	if(mexp==0) 		/* subnormal number */
	  {pt[n0]=B0;		/* scale t to be normal, 2^54 or 2^18 */
	   t*=x; pt[n0]=pt[n0]/3+B2;
	  }
	else
	  pt[n0]=px[n0]/3+B1;

    /* new cbrt to 23 bits, may be implemented in single precision */
	r=t*t/x;
	s=C+r*t;
	t*=G+F/(s+E+D/s);

    /* chopped to 20 bits and make it larger than cbrt(x) */
	pt[n1]=0; pt[n0]+=0x00000001;

    /* one step newton iteration to 53 bits with error less than 0.667 ulps */
	s=t*t;		/* t*t is exact */
	r=x/s;
	w=t+t;
	r=(r-t)/(w+r);	/* r-s is exact */
	t=t+t*r;

    /* restore the sign bit */
	pt[n0] |= sign;
	return(t);
    }

/*******************************************************************************
*
* mathSoftCosh - software emulation for floating point hyperbolic cosine
*
*
*	Computes hyperbolic cosine from:
*
*		cosh(X) = 0.5 * (exp(X) + exp(-X))
*
*	Returns double precision hyperbolic cosine of double precision
*	floating point number.
*
*	Inputs greater than log(DBL_MAX) result in overflow.
*	Inputs less than log(DBL_MIN) result in underflow.
*
*	For precision information refer to documentation of the
*	floating point library routines called.
*
* RETURNS: double-precision hyperbolic cosine
*
* AUTHOR: Fred Fish
*
*				N O T I C E
*
*			Copyright Abandoned, 1987, Fred Fish
*
*	This previously copyrighted work has been placed into the
*	public domain by the author (Fred Fish) and may be freely used
*	for any purpose, private or commercial.  I would appreciate
*	it, as a courtesy, if this notice is left in all copies and
*	derivative works.  Thank you, and enjoy...
*
*/

LOCAL double mathSoftCosh 
    (
    double	dblParam 
    )
    {
    double	retValue;

    if (dblParam > LOGE_MAXDOUBLE)
	{
	retValue = HUGE_VAL;
        }
    else if (dblParam < LOGE_MINDOUBLE)
	{
	retValue = -HUGE_VAL;
    	}
    else
	{
	dblParam = exp (dblParam);
	retValue = 0.5 * (dblParam + 1.0/dblParam);
    	}

    return (retValue);
    }

/*******************************************************************************
*
* mathSoftSinh - software emulation for floating point hyperbolic sine
*
*	Computes hyperbolic sine from:
*
*		sinh(x) = 0.5 * (exp(x) - 1.0/exp(x))
*
*
*	Returns double precision hyperbolic sine of double precision
*	floating point number.
*
*	Inputs greater than ln(DBL_MAX) result in overflow.
*	Inputs less than ln(DBL_MIN) result in underflow.
*
* RETURNS: double-precision hyperbolic sine
*
* AUTHOR: Fred Fish
*
*				N O T I C E
*
*			Copyright Abandoned, 1987, Fred Fish
*
*	This previously copyrighted work has been placed into the
*	public domain by the author (Fred Fish) and may be freely used
*	for any purpose, private or commercial.  I would appreciate
*	it, as a courtesy, if this notice is left in all copies and
*	derivative works.  Thank you, and enjoy...
*
*/

LOCAL double mathSoftSinh 
    (
    double	dblParam 
    )
    {
    double	retValue;


    if (dblParam > LOGE_MAXDOUBLE)
	{
	retValue = HUGE_VAL;
    	}
    else if (dblParam < LOGE_MINDOUBLE)
	{
	retValue = -HUGE_VAL;
    	}
    else
	{
	dblParam = exp (dblParam);
	retValue = 0.5 * (dblParam - (1.0 / dblParam));
    	}

    return (retValue);
    }

/*******************************************************************************
*
* mathSoftTanh - software emulation for floating point hyperbolic tangent
*
*
*	Computes hyperbolic tangent from:
*
*		tanh(x) = 1.0 for x > TANH_MAXARG
*			= -1.0 for x < -TANH_MAXARG
*			=  sinh(x) / cosh(x) otherwise
*
*	Returns double precision hyperbolic tangent of double precision
*	floating point number.
*
* RETURNS: double-precision hyperbolic tangent
*
* AUTHOR: Fred Fish
*
*				N O T I C E
*
*			Copyright Abandoned, 1987, Fred Fish
*
*	This previously copyrighted work has been placed into the
*	public domain by the author (Fred Fish) and may be freely used
*	for any purpose, private or commercial.  I would appreciate
*	it, as a courtesy, if this notice is left in all copies and
*	derivative works.  Thank you, and enjoy...
*/

LOCAL double mathSoftTanh 
    (
    double	dblParam 
    )
    {
    double	retValue;

    if (dblParam > TANH_MAXARG || dblParam < -TANH_MAXARG)
	{
	if (dblParam > 0.0)
	    retValue = 1.0;
	else
	    retValue = -1.0;
    	}
    else
	{
	retValue = sinh (dblParam) / cosh (dblParam);
    	}

    return (retValue);
    }

/*******************************************************************************
*
* mathSoftLog2 - software emulation for floating point log base 2
*
* This routine returns a double-precision log base 2 of a double-precision
* floating point number.
*
*	Computes log2(x) from:
*
*		log2(x) = log2(e) * log(x)
*
*
* RETURNS: double-precision log base 2
*
*/

LOCAL double mathSoftLog2 
    (
    double	dblParam 
    )
    {

    return (LOG2E * log (dblParam));
    }

/*******************************************************************************
*
* mathSoftSincos - software emulation for simultaneous sine and cosine
*
* This routine obtains both the sine and cosine for the specified
* floating point value and returns both.
*
* NOMANUAL
*/

void mathSoftSincos 
    (
    double	dblParam,		/* angle in radians */
    double	*pSinResult,		/* sine result buffer */
    double	*pCosResult 		/* cosine result buffer */
    )
    {

    *pSinResult = sin (dblParam);
    *pCosResult = cos (dblParam);
    }


/******************************************************************************
*
* mathSoftFabsf - single-precision software floating point absolute value
*
* This routine takes the input single-precision floating point
* parameter and returns the absolute value.
*
* RETURNS: single-precision absolute value.
*/

LOCAL float mathSoftFabsf 
    (
    float	fltParam 
    )
    {

    return ((fltParam < 0.0) ? -fltParam : fltParam);
    }

/******************************************************************************
*
* mathSoftCeilf - single-precision software floating point ceiling
*
* This routine takes the input single-precision floating point
* parameter and returns (in single-precision floating point format)
* the integer immediately greater than the input parameter.
*
* RETURNS: single-precision representation of next largest integer.
*/

LOCAL float mathSoftCeilf 
    (
    float	fltParam 
    )
    {

    return (-floorf (-fltParam));
    }

/******************************************************************************
*
* mathSoftPowf - single-precision software floating point power function
*
* This routine takes two input single-precision floating point
* parameters, <fltX> and <fltY>, and returns the single-precision
* value of <fltX> to the <fltY> power.
*
* INTERNAL
* The US Software emulation library has a special function for taking
* a floating point number to an integer power.  This routine therefore
* checks to see if the <fltY> parameter is an integer, and uses the
* US Software function (XTOI) if it is.
*
* RETURNS: single-precision value of <fltX> to <fltY> power.
*/

LOCAL float mathSoftPowf 
    (
    float	fltX,
    float	fltY 
    )
    {

    if (isNaN_SINGLE (fltY))
	return (fltY);			/* fltY = NaN --> NaN */

    if (fltX == 1.0)
	return (1.0);			/* fltX = 1 --> 1 */

    if (fltY == floorf (fltY))		/* int fltY --> XTOI(fltX,fltY) */
	return (mathSoftRealtointf (fltX, (long int) fltY));

    return (expf (fltY * logf (fltX)));
    }