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)));    }