OpenCores

/w_cabs.c

0,0 → 1,20

/sf_isinf.c

0,0 → 1,48

/e_acosh.c

0,0 → 1,135

/Makefile.in

0,0 → 1,560

/s_mathcnst.c

0,0 → 1,24

/s_floor.c

0,0 → 1,92

/sf_erf.c

0,0 → 1,246

/wf_cabs.c

0,0 → 1,20

/s_infconst.c

0,0 → 1,15

/s_fmod.c

0,0 → 1,187

/ef_acosh.c

0,0 → 1,53

/e_j0.c

0,0 → 1,487

/sf_fmod.c

0,0 → 1,103

/s_exp.c

0,0 → 1,133

/erf_gamma.c

0,0 → 1,34

/s_sine.c

0,0 → 1,166

/s_numtest.c

0,0 → 1,58

/ef_j1.c

0,0 → 1,439

/sf_sine.c

0,0 → 1,112

sf_sine.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: s_atan2.c =================================================================== --- s_atan2.c (nonexistent) +++ s_atan2.c (revision 1765) @@ -0,0 +1,89 @@ + +/* @(#)z_atan2.c 1.0 98/08/13 */ + +/* +FUNCTION + <>, <>---arc tangent of y/x + +INDEX + atan2 +INDEX + atan2f + +ANSI_SYNOPSIS + #include + double atan2(double <[y]>,double <[x]>); + float atan2f(float <[y]>,float <[x]>); + +TRAD_SYNOPSIS + #include + double atan2(<[y]>,<[x]>); + double <[y]>; + double <[x]>; + + float atan2f(<[y]>,<[x]>); + float <[y]>; + float <[x]>; + +DESCRIPTION + +<> computes the inverse tangent (arc tangent) of <[y]>/<[x]>. +<> produces the correct result even for angles near +@ifinfo +pi/2 or -pi/2 +@end ifinfo +@tex +$\pi/2$ or $-\pi/2$ +@end tex +(that is, when <[x]> is near 0). + +<> is identical to <>, save that it takes and returns +<>. + +RETURNS +<> and <> return a value in radians, in the range of +@ifinfo +-pi to pi. +@end ifinfo +@tex +$-\pi$ to $\pi$. +@end tex + +If both <[x]> and <[y]> are 0.0, <> causes a <> error. + +You can modify error handling for these functions using <>. + +PORTABILITY +<> is ANSI C. <> is an extension. + + +*/ + +/****************************************************************** + * Arctangent2 + * + * Input: + * v, u - floating point values + * + * Output: + * arctan2 of v / u + * + * Description: + * This routine returns the arctan2 of v / u. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +double +_DEFUN (atan2, (double, double), + double v _AND + double u) +{ + return (atangent (0.0, v, u, 1)); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: s_sineh.c =================================================================== --- s_sineh.c (nonexistent) +++ s_sineh.c (revision 1765) @@ -0,0 +1,185 @@ + +/* @(#)z_sineh.c 1.0 98/08/13 */ +/****************************************************************** + * The following routines are coded directly from the algorithms + * and coefficients given in "Software Manual for the Elementary + * Functions" by William J. Cody, Jr. and William Waite, Prentice + * Hall, 1980. + ******************************************************************/ + +/* +FUNCTION + <>, <>, <>, <>, <>---hyperbolic sine or cosine + +INDEX + sinh +INDEX + sinhf +INDEX + cosh +INDEX + coshf + +ANSI_SYNOPSIS + #include + double sinh(double <[x]>); + float sinhf(float <[x]>); + double cosh(double <[x]>); + float coshf(float <[x]>); +TRAD_SYNOPSIS + #include + double sinh(<[x]>) + double <[x]>; + + float sinhf(<[x]>) + float <[x]>; + + double cosh(<[x]>) + double <[x]>; + + float coshf(<[x]>) + float <[x]>; + +DESCRIPTION + <> and <> compute the hyperbolic sine or cosine + of the argument <[x]>. + Angles are specified in radians. <>(<[x]>) is defined as + @ifinfo + . (exp(<[x]>) - exp(-<[x]>))/2 + @end ifinfo + @tex + $${e^x - e^{-x}}\over 2$$ + @end tex + <> is defined as + @ifinfo + . (exp(<[x]>) - exp(-<[x]>))/2 + @end ifinfo + @tex + $${e^x + e^{-x}}\over 2$$ + @end tex + + <> and <> are identical, save that they take + and returns <> values. + +RETURNS + The hyperbolic sine or cosine of <[x]> is returned. + + When the correct result is too large to be representable (an + overflow), the functions return <> with the + appropriate sign, and sets the global value <> to + <>. + +PORTABILITY + <> is ANSI C. + <> is an extension. + <> is ANSI C. + <> is an extension. + +*/ + +/****************************************************************** + * Hyperbolic Sine + * + * Input: + * x - floating point value + * + * Output: + * hyperbolic sine of x + * + * Description: + * This routine calculates hyperbolic sines. + * + *****************************************************************/ + +#include +#include "fdlibm.h" +#include "zmath.h" + +static const double q[] = { -0.21108770058106271242e+7, + 0.36162723109421836460e+5, + -0.27773523119650701667e+3 }; +static const double p[] = { -0.35181283430177117881e+6, + -0.11563521196851768270e+5, + -0.16375798202630751372e+3, + -0.78966127417357099479 }; +static const double LNV = 0.6931610107421875000; +static const double INV_V2 = 0.24999308500451499336; +static const double V_OVER2_MINUS1 = 0.13830277879601902638e-4; + +double +_DEFUN (sineh, (double, int), + double x _AND + int cosineh) +{ + double y, f, P, Q, R, res, z, w; + int sgn = 1; + double WBAR = 18.55; + + /* Check for special values. */ + switch (numtest (x)) + { + case NAN: + errno = EDOM; + return (x); + case INF: + errno = ERANGE; + return (ispos (x) ? z_infinity.d : -z_infinity.d); + } + + y = fabs (x); + + if (!cosineh && x < 0.0) + sgn = -1; + + if ((y > 1.0 && !cosineh) || cosineh) + { + if (y > BIGX) + { + w = y - LNV; + + /* Check for w > maximum here. */ + if (w > BIGX) + { + errno = ERANGE; + return (x); + } + + z = exp (w); + + if (w > WBAR) + res = z * (V_OVER2_MINUS1 + 1.0); + } + + else + { + z = exp (y); + if (cosineh) + res = (z + 1 / z) / 2.0; + else + res = (z - 1 / z) / 2.0; + } + + if (sgn < 0) + res = -res; + } + else + { + /* Check for y being too small. */ + if (y < z_rooteps) + { + res = x; + } + /* Calculate the Taylor series. */ + else + { + f = x * x; + Q = ((f + q[2]) * f + q[1]) * f + q[0]; + P = ((p[3] * f + p[2]) * f + p[1]) * f + p[0]; + R = f * (P / Q); + + res = x + x * R; + } + } + + return (res); +} Index: zmath.h =================================================================== --- zmath.h (nonexistent) +++ zmath.h (revision 1765) @@ -0,0 +1,55 @@ +#ifndef __ZMATH_H +#define __ZMATH_H + +#include + +#define NUM 3 +#define NAN 2 +#define INF 1 + +#define __PI 3.14159265358979323846 +#define __SQRT_HALF 0.70710678118654752440 +#define __PI_OVER_TWO 1.57079632679489661923132 + +extern double BIGX; +extern double SMALLX; + +typedef const union +{ + long l[2]; + double d; +} udouble; + +typedef const union +{ + long l; + float f; +} ufloat; + +extern double BIGX; +extern double SMALLX; + +extern udouble z_infinity; +extern udouble z_notanum; +extern double z_rooteps; + +extern ufloat z_infinity_f; +extern ufloat z_notanum_f; +extern float z_rooteps_f; + +/* Core math routines. */ + +int _EXFUN (numtest, (double)); +int _EXFUN (numtestf, (float)); +double _EXFUN (logarithm, (double, int)); +float _EXFUN (logarithmf, (float, int)); +double _EXFUN (sine, (double, int)); +float _EXFUN (sinef, (float, int)); +double _EXFUN (asine, (double, int)); +float _EXFUN (asinef, (float, int)); +double _EXFUN (atangent, (double, double, double, int)); +float _EXFUN (atangentf, (float, float, float, int)); +double _EXFUN (sineh, (double, int)); +float _EXFUN (sinehf, (float, int)); + +#endif /* no __ZMATH_H */ Index: w_jn.c =================================================================== --- w_jn.c (nonexistent) +++ w_jn.c (revision 1765) @@ -0,0 +1,248 @@ + +/* @(#)w_jn.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* +FUNCTION +<>,<>,<>,<>---Bessel functions + +INDEX +j0 +INDEX +j0f +INDEX +j1 +INDEX +j1f +INDEX +jn +INDEX +jnf +INDEX +y0 +INDEX +y0f +INDEX +y1 +INDEX +y1f +INDEX +yn +INDEX +ynf + +ANSI_SYNOPSIS +#include +double j0(double <[x]>); +float j0f(float <[x]>); +double j1(double <[x]>); +float j1f(float <[x]>); +double jn(int <[n]>, double <[x]>); +float jnf(int <[n]>, float <[x]>); +double y0(double <[x]>); +float y0f(float <[x]>); +double y1(double <[x]>); +float y1f(float <[x]>); +double yn(int <[n]>, double <[x]>); +float ynf(int <[n]>, float <[x]>); + +TRAD_SYNOPSIS +#include + +double j0(<[x]>) +double <[x]>; +float j0f(<[x]>) +float <[x]>; +double j1(<[x]>) +double <[x]>; +float j1f(<[x]>) +float <[x]>; +double jn(<[n]>, <[x]>) +int <[n]>; +double <[x]>; +float jnf(<[n]>, <[x]>) +int <[n]>; +float <[x]>; + +double y0(<[x]>) +double <[x]>; +float y0f(<[x]>) +float <[x]>; +double y1(<[x]>) +double <[x]>; +float y1f(<[x]>) +float <[x]>; +double yn(<[n]>, <[x]>) +int <[n]>; +double <[x]>; +float ynf(<[n]>, <[x]>) +int <[n]>; +float <[x]>; + +DESCRIPTION +The Bessel functions are a family of functions that solve the +differential equation +@ifinfo +. 2 2 2 +. x y'' + xy' + (x - p )y = 0 +@end ifinfo +@tex +$$x^2{d^2y\over dx^2} + x{dy\over dx} + (x^2-p^2)y = 0$$ +@end tex +These functions have many applications in engineering and physics. + +<> calculates the Bessel function of the first kind of order +<[n]>. <> and <> are special cases for order 0 and order +1 respectively. + +Similarly, <> calculates the Bessel function of the second kind of +order <[n]>, and <> and <> are special cases for order 0 and +1. + +<>, <>, <>, <>, <>, and <> perform the +same calculations, but on <> rather than <> values. + +RETURNS +The value of each Bessel function at <[x]> is returned. + +PORTABILITY +None of the Bessel functions are in ANSI C. +*/ + +/* + * wrapper jn(int n, double x), yn(int n, double x) + * floating point Bessel's function of the 1st and 2nd kind + * of order n + * + * Special cases: + * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal; + * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal. + * Note 2. About jn(n,x), yn(n,x) + * For n=0, j0(x) is called, + * for n=1, j1(x) is called, + * for nx, a continued fraction approximation to + * j(n,x)/j(n-1,x) is evaluated and then backward + * recursion is used starting from a supposed value + * for j(n,x). The resulting value of j(0,x) is + * compared with the actual value to correct the + * supposed value of j(n,x). + * + * yn(n,x) is similar in all respects, except + * that forward recursion is used for all + * values of n>1. + * + */ + +#include "fdlibm.h" +#include + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double jn(int n, double x) /* wrapper jn */ +#else + double jn(n,x) /* wrapper jn */ + double x; int n; +#endif +{ +#ifdef _IEEE_LIBM + return jn(n,x); +#else + double z; + struct exception exc; + z = jn(n,x); + if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z; + if(fabs(x)>X_TLOSS) { + /* jn(|x|>X_TLOSS) */ + exc.type = TLOSS; + exc.name = "jn"; + exc.err = 0; + exc.arg1 = n; + exc.arg2 = x; + exc.retval = 0.0; + if (_LIB_VERSION == _POSIX_) + errno = ERANGE; + else if (!matherr(&exc)) { + errno = ERANGE; + } + if (exc.err != 0) + errno = exc.err; + return exc.retval; + } else + return z; +#endif +} + +#ifdef __STDC__ + double yn(int n, double x) /* wrapper yn */ +#else + double yn(n,x) /* wrapper yn */ + double x; int n; +#endif +{ +#ifdef _IEEE_LIBM + return yn(n,x); +#else + double z; + struct exception exc; + z = yn(n,x); + if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z; + if(x <= 0.0){ + /* yn(n,0) = -inf or yn(x<0) = NaN */ +#ifndef HUGE_VAL +#define HUGE_VAL inf + double inf = 0.0; + + SET_HIGH_WORD(inf,0x7ff00000); /* set inf to infinite */ +#endif + exc.type = DOMAIN; /* should be SING for IEEE */ + exc.name = "yn"; + exc.err = 0; + exc.arg1 = n; + exc.arg2 = x; + if (_LIB_VERSION == _SVID_) + exc.retval = -HUGE; + else + exc.retval = -HUGE_VAL; + if (_LIB_VERSION == _POSIX_) + errno = EDOM; + else if (!matherr(&exc)) { + errno = EDOM; + } + if (exc.err != 0) + errno = exc.err; + return exc.retval; + } + if(x>X_TLOSS) { + /* yn(x>X_TLOSS) */ + exc.type = TLOSS; + exc.name = "yn"; + exc.err = 0; + exc.arg1 = n; + exc.arg2 = x; + exc.retval = 0.0; + if (_LIB_VERSION == _POSIX_) + errno = ERANGE; + else if (!matherr(&exc)) { + errno = ERANGE; + } + if (exc.err != 0) + errno = exc.err; + return exc.retval; + } else + return z; +#endif +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ Index: sf_isnan.c =================================================================== --- sf_isnan.c (nonexistent) +++ sf_isnan.c (revision 1765) @@ -0,0 +1,48 @@ + +/* @(#)z_isnanf.c 1.0 98/08/13 */ +/****************************************************************** + * isnanf + * + * Input: + * x - pointer to a floating point value + * + * Output: + * An integer that indicates if the number is NaN. + * + * Description: + * This routine returns an integer that indicates if the number + * passed in is NaN (1) or is finite (0). + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +int +_DEFUN (isnanf, (float), + float x) +{ + __int32_t wx; + int exp; + + GET_FLOAT_WORD (wx, x); + exp = (wx & 0x7f800000) >> 23; + + if ((exp == 0x7f8) && (wx & 0x7fffff)) + return (1); + else + return (0); +} + + +#ifdef _DOUBLE_IS_32BITS + +int +_DEFUN (isnan, (double), + double x) +{ + return isnanf((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ + Index: s_signif.c =================================================================== --- s_signif.c (nonexistent) +++ s_signif.c (revision 1765) @@ -0,0 +1,34 @@ + +/* @(#)s_signif.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * significand(x) computes just + * scalb(x, (double) -ilogb(x)), + * for exercising the fraction-part(F) IEEE 754-1985 test vector. + */ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double significand(double x) +#else + double significand(x) + double x; +#endif +{ + return scalb(x,(double) -ilogb(x)); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: wf_jn.c =================================================================== --- wf_jn.c (nonexistent) +++ wf_jn.c (revision 1765) @@ -0,0 +1,138 @@ +/* wf_jn.c -- float version of w_jn.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" +#include + + +#ifdef __STDC__ + float jnf(int n, float x) /* wrapper jnf */ +#else + float jnf(n,x) /* wrapper jnf */ + float x; int n; +#endif +{ +#ifdef _IEEE_LIBM + return jnf(n,x); +#else + float z; + struct exception exc; + z = jnf(n,x); + if(_LIB_VERSION == _IEEE_ || isnanf(x) ) return z; + if(fabsf(x)>(float)X_TLOSS) { + /* jnf(|x|>X_TLOSS) */ + exc.type = TLOSS; + exc.name = "jnf"; + exc.err = 0; + exc.arg1 = (double)n; + exc.arg2 = (double)x; + exc.retval = 0.0; + if (_LIB_VERSION == _POSIX_) + errno = ERANGE; + else if (!matherr(&exc)) { + errno = ERANGE; + } + if (exc.err != 0) + errno = exc.err; + return exc.retval; + } else + return z; +#endif +} + +#ifdef __STDC__ + float ynf(int n, float x) /* wrapper ynf */ +#else + float ynf(n,x) /* wrapper ynf */ + float x; int n; +#endif +{ +#ifdef _IEEE_LIBM + return ynf(n,x); +#else + float z; + struct exception exc; + z = ynf(n,x); + if(_LIB_VERSION == _IEEE_ || isnanf(x) ) return z; + if(x <= (float)0.0){ + /* ynf(n,0) = -inf or ynf(x<0) = NaN */ +#ifndef HUGE_VAL +#define HUGE_VAL inf + double inf = 0.0; + + SET_HIGH_WORD(inf,0x7ff00000); /* set inf to infinite */ +#endif + exc.type = DOMAIN; /* should be SING for IEEE */ + exc.name = "ynf"; + exc.err = 0; + exc.arg1 = (double)n; + exc.arg2 = (double)x; + if (_LIB_VERSION == _SVID_) + exc.retval = -HUGE; + else + exc.retval = -HUGE_VAL; + if (_LIB_VERSION == _POSIX_) + errno = EDOM; + else if (!matherr(&exc)) { + errno = EDOM; + } + if (exc.err != 0) + errno = exc.err; + return (float)exc.retval; + } + if(x>(float)X_TLOSS) { + /* ynf(x>X_TLOSS) */ + exc.type = TLOSS; + exc.name = "ynf"; + exc.err = 0; + exc.arg1 = (double)n; + exc.arg2 = (double)x; + exc.retval = 0.0; + if (_LIB_VERSION == _POSIX_) + errno = ERANGE; + else if (!matherr(&exc)) { + errno = ERANGE; + } + if (exc.err != 0) + errno = exc.err; + return (float)exc.retval; + } else + return z; +#endif +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double jn(int n, double x) +#else + double jn(n,x) + double x; int n; +#endif +{ + return (double) jnf(n, (float) x); +} + +#ifdef __STDC__ + double yn(int n, double x) +#else + double yn(n,x) + double x; int n; +#endif +{ + return (double) ynf(n, (float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ Index: s_log10.c =================================================================== --- s_log10.c (nonexistent) +++ s_log10.c (revision 1765) @@ -0,0 +1,68 @@ + +/* @(#)z_log10.c 1.0 98/08/13 */ +/****************************************************************** + * Logarithm + * + * Input: + * x - floating point value + * + * Output: + * logarithm of x + * + * Description: + * This routine returns the logarithm of x (base 10). + * + *****************************************************************/ + +/* +FUNCTION + <>, <>---base 10 logarithms + +INDEX +log10 +INDEX +log10f + +ANSI_SYNOPSIS + #include + double log10(double <[x]>); + float log10f(float <[x]>); + +TRAD_SYNOPSIS + #include + double log10(<[x]>) + double <[x]>; + + float log10f(<[x]>) + float <[x]>; + +DESCRIPTION +<> returns the base 10 logarithm of <[x]>. +It is implemented as <) / log(10)>>. + +<> is identical, save that it takes and returns <> values. + +RETURNS +<> and <> return the calculated value. + +See the description of <> for information on errors. + +PORTABILITY +<> is ANSI C. <> is an extension. + +*/ + + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +double +_DEFUN (log10, (double), + double x) +{ + return (logarithm (x, 1)); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: sf_ldexp.c =================================================================== --- sf_ldexp.c (nonexistent) +++ sf_ldexp.c (revision 1765) @@ -0,0 +1,81 @@ + +/* @(#)z_ldexpf.c 1.0 98/08/13 */ +/****************************************************************** + * ldexp + * + * Input: + * d - a floating point value + * e - an exponent value + * + * Output: + * A floating point value f such that f = d * 2 ^ e. + * + * Description: + * This function creates a floating point number f such that + * f = d * 2 ^ e. + * + *****************************************************************/ + +#include +#include "fdlibm.h" +#include "zmath.h" + +#define FLOAT_EXP_OFFS 127 + +float +_DEFUN (ldexpf, (float, int), + float d _AND + int e) +{ + int exp; + __int32_t wd; + + GET_FLOAT_WORD (wd, d); + + /* Check for special values and then scale d by e. */ + switch (numtestf (wd)) + { + case NAN: + errno = EDOM; + break; + + case INF: + errno = ERANGE; + break; + + case 0: + break; + + default: + exp = (wd & 0x7f800000) >> 23; + exp += e; + + if (exp > FLT_MAX_EXP + FLOAT_EXP_OFFS) + { + errno = ERANGE; + d = z_infinity_f.f; + } + else if (exp < FLT_MIN_EXP + FLOAT_EXP_OFFS) + { + errno = ERANGE; + d = -z_infinity_f.f; + } + else + { + wd &= 0x807fffff; + wd |= exp << 23; + SET_FLOAT_WORD (d, wd); + } + } + + return (d); +} + +#ifdef _DOUBLE_IS_32BITS + +double ldexp (double x, int e) +{ + return (double) ldexpf ((float) x, e); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ Index: s_atan.c =================================================================== --- s_atan.c (nonexistent) +++ s_atan.c (revision 1765) @@ -0,0 +1,83 @@ + +/* @(#)z_atan.c 1.0 98/08/13 */ + +/* +FUNCTION + <>, <>---arc tangent + +INDEX + atan +INDEX + atanf + +ANSI_SYNOPSIS + #include + double atan(double <[x]>); + float atanf(float <[x]>); + +TRAD_SYNOPSIS + #include + double atan(<[x]>); + double <[x]>; + + float atanf(<[x]>); + float <[x]>; + +DESCRIPTION + +<> computes the inverse tangent (arc tangent) of the input value. + +<> is identical to <>, save that it operates on <>. + +RETURNS +@ifinfo +<> returns a value in radians, in the range of -pi/2 to pi/2. +@end ifinfo +@tex +<> returns a value in radians, in the range of $-\pi/2$ to $\pi/2$. +@end tex + +PORTABILITY +<> is ANSI C. <> is an extension. + +*/ + +/****************************************************************** + * Arctangent + * + * Input: + * x - floating point value + * + * Output: + * arctan of x + * + * Description: + * This routine returns the arctan of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +double +_DEFUN (atan, (double), + double x) +{ + switch (numtest (x)) + { + case NAN: + errno = EDOM; + return (x); + case INF: + /* this should check to see if neg NaN or pos NaN... */ + return (__PI_OVER_TWO); + case 0: + return (0.0); + default: + return (atangent (x, 0, 0, 0)); + } +} + +#endif /* _DOUBLE_IS_32BITS */ Index: s_frexp.c =================================================================== --- s_frexp.c (nonexistent) +++ s_frexp.c (revision 1765) @@ -0,0 +1,110 @@ + +/* @(#)z_frexp.c 1.0 98/08/13 */ + +/* +FUNCTION + <>, <>---split floating-point number +INDEX + frexp +INDEX + frexpf + +ANSI_SYNOPSIS + #include + double frexp(double <[val]>, int *<[exp]>); + float frexpf(float <[val]>, int *<[exp]>); + +TRAD_SYNOPSIS + #include + double frexp(<[val]>, <[exp]>) + double <[val]>; + int *<[exp]>; + + float frexpf(<[val]>, <[exp]>) + float <[val]>; + int *<[exp]>; + + +DESCRIPTION + All non zero, normal numbers can be described as <[m]> * 2**<[p]>. + <> represents the double <[val]> as a mantissa <[m]> + and a power of two <[p]>. The resulting mantissa will always + be greater than or equal to <<0.5>>, and less than <<1.0>> (as + long as <[val]> is nonzero). The power of two will be stored + in <<*>><[exp]>. + +@ifinfo +<[m]> and <[p]> are calculated so that +<[val]> is <[m]> times <<2>> to the power <[p]>. +@end ifinfo +@tex +<[m]> and <[p]> are calculated so that +$ val = m \times 2^p $. +@end tex + +<> is identical, other than taking and returning +floats rather than doubles. + +RETURNS +<> returns the mantissa <[m]>. If <[val]> is <<0>>, infinity, +or Nan, <> will set <<*>><[exp]> to <<0>> and return <[val]>. + +PORTABILITY +<> is ANSI. +<> is an extension. + + +*/ + +/***************************************************************** + * frexp + * + * Input: + * d - floating point value + * exp - exponent value + * + * Output: + * A floating point value in the range [0.5, 1). + * + * Description: + * This routine breaks a floating point value into a number f and + * an exponent exp such that d = f * 2 ^ exp. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +double frexp (double d, int *exp) +{ + double f; + __uint32_t hd, ld, hf, lf; + + EXTRACT_WORDS (hd, ld, d); + + /* Get the exponent. */ + *exp = ((hd & 0x7ff00000) >> 20) - 1022; + + /* Get the mantissa. */ + lf = ld; + hf = hd & 0x800fffff; + hf |= 0x3fe00000; + + INSERT_WORDS (f, hf, lf); + + /* Check for special values. */ + switch (numtest (f)) + { + case NAN: + case INF: + errno = EDOM; + *exp = 0; + return (f); + } + + return (f); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: s_acos.c =================================================================== --- s_acos.c (nonexistent) +++ s_acos.c (revision 1765) @@ -0,0 +1,93 @@ + +/* @(#)z_acos.c 1.0 98/08/13 */ + +/* +FUNCTION + <>, <>---arc cosine + +INDEX + acos +INDEX + acosf + +ANSI_SYNOPSIS + #include + double acos(double <[x]>); + float acosf(float <[x]>); + +TRAD_SYNOPSIS + #include + double acos(<[x]>) + double <[x]>; + + float acosf(<[x]>) + float <[x]>; + + + +DESCRIPTION + + <> computes the inverse cosine (arc cosine) of the input value. + Arguments to <> must be in the range @minus{}1 to 1. + + <> is identical to <>, except that it performs + its calculations on <>. + +RETURNS + @ifinfo + <> and <> return values in radians, in the range of 0 to pi +. + @end ifinfo + @tex + <> and <> return values in radians, in the range of <<0>> t +o $\pi$. + @end tex + + If <[x]> is not between @minus{}1 and 1, the returned value is NaN + (not a number) the global variable <> is set to <>, and a + <> message is sent as standard error output. + + You can modify error handling for these functions using <>. + + +QUICKREF ANSI SVID POSIX RENTRANT + acos y,y,y,m + acosf n,n,n,m + +MATHREF + acos, [-1,1], acos(arg),,, + acos, NAN, arg,DOMAIN,EDOM + +MATHREF + acosf, [-1,1], acosf(arg),,, + acosf, NAN, argf,DOMAIN,EDOM + +*/ + +/***************************************************************** + * Arccosine + * + * Input: + * x - floating point value + * + * Output: + * arccosine of x + * + * Description: + * This routine returns the arccosine of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +double +_DEFUN (acos, (double), + double x) +{ + return (asine (x, 1)); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: sf_signif.c =================================================================== --- sf_signif.c (nonexistent) +++ sf_signif.c (revision 1765) @@ -0,0 +1,40 @@ +/* sf_signif.c -- float version of s_signif.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ + float significandf(float x) +#else + float significandf(x) + float x; +#endif +{ + return scalbf(x,(float) -ilogbf(x)); +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double significand(double x) +#else + double significand(x) + double x; +#endif +{ + return (double) significandf((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ Index: e_hypot.c =================================================================== --- e_hypot.c (nonexistent) +++ e_hypot.c (revision 1765) @@ -0,0 +1,170 @@ + +/* @(#)e_hypot.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* +FUNCTION + <>, <>---distance from origin +INDEX + hypot +INDEX + hypotf + +ANSI_SYNOPSIS + #include + double hypot(double <[x]>, double <[y]>); + float hypotf(float <[x]>, float <[y]>); + +TRAD_SYNOPSIS + double hypot(<[x]>, <[y]>) + double <[x]>, <[y]>; + + float hypotf(<[x]>, <[y]>) + float <[x]>, <[y]>; + +DESCRIPTION + <> calculates the Euclidean distance + @tex + $\sqrt{x^2+y^2}$ + @end tex + @ifinfo + <*<[x]> + <[y]>*<[y]>)>> + @end ifinfo + between the origin (0,0) and a point represented by the + Cartesian coordinates (<[x]>,<[y]>). <> differs only + in the type of its arguments and result. + +RETURNS + Normally, the distance value is returned. On overflow, + <> returns <> and sets <> to + <>. + + You can change the error treatment with <>. + +PORTABILITY + <> and <> are not ANSI C. */ + +/* hypot(x,y) + * + * Method : + * If (assume round-to-nearest) z=x*x+y*y + * has error less than sqrt(2)/2 ulp, than + * sqrt(z) has error less than 1 ulp (exercise). + * + * So, compute sqrt(x*x+y*y) with some care as + * follows to get the error below 1 ulp: + * + * Assume x>y>0; + * (if possible, set rounding to round-to-nearest) + * 1. if x > 2y use + * x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y + * where x1 = x with lower 32 bits cleared, x2 = x-x1; else + * 2. if x <= 2y use + * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y)) + * where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1, + * y1= y with lower 32 bits chopped, y2 = y-y1. + * + * NOTE: scaling may be necessary if some argument is too + * large or too tiny + * + * Special cases: + * hypot(x,y) is INF if x or y is +INF or -INF; else + * hypot(x,y) is NAN if x or y is NAN. + * + * Accuracy: + * hypot(x,y) returns sqrt(x^2+y^2) with error less + * than 1 ulps (units in the last place) + */ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double hypot(double x, double y) +#else + double hypot(x,y) + double x, y; +#endif +{ + double a=x,b=y,t1,t2,y1,y2,w; + __int32_t j,k,ha,hb; + + GET_HIGH_WORD(ha,x); + ha &= 0x7fffffff; + GET_HIGH_WORD(hb,y); + hb &= 0x7fffffff; + if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;} + SET_HIGH_WORD(a,ha); /* a <- |a| */ + SET_HIGH_WORD(b,hb); /* b <- |b| */ + if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */ + k=0; + if(ha > 0x5f300000) { /* a>2**500 */ + if(ha >= 0x7ff00000) { /* Inf or NaN */ + __uint32_t low; + w = a+b; /* for sNaN */ + GET_LOW_WORD(low,a); + if(((ha&0xfffff)|low)==0) w = a; + GET_LOW_WORD(low,b); + if(((hb^0x7ff00000)|low)==0) w = b; + return w; + } + /* scale a and b by 2**-600 */ + ha -= 0x25800000; hb -= 0x25800000; k += 600; + SET_HIGH_WORD(a,ha); + SET_HIGH_WORD(b,hb); + } + if(hb < 0x20b00000) { /* b < 2**-500 */ + if(hb <= 0x000fffff) { /* subnormal b or 0 */ + __uint32_t low; + GET_LOW_WORD(low,b); + if((hb|low)==0) return a; + t1=0; + SET_HIGH_WORD(t1,0x7fd00000); /* t1=2^1022 */ + b *= t1; + a *= t1; + k -= 1022; + } else { /* scale a and b by 2^600 */ + ha += 0x25800000; /* a *= 2^600 */ + hb += 0x25800000; /* b *= 2^600 */ + k -= 600; + SET_HIGH_WORD(a,ha); + SET_HIGH_WORD(b,hb); + } + } + /* medium size a and b */ + w = a-b; + if (w>b) { + t1 = 0; + SET_HIGH_WORD(t1,ha); + t2 = a-t1; + w = sqrt(t1*t1-(b*(-b)-t2*(a+t1))); + } else { + a = a+a; + y1 = 0; + SET_HIGH_WORD(y1,hb); + y2 = b - y1; + t1 = 0; + SET_HIGH_WORD(t1,ha+0x00100000); + t2 = a - t1; + w = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b))); + } + if(k!=0) { + __uint32_t high; + t1 = 1.0; + GET_HIGH_WORD(high,t1); + SET_HIGH_WORD(t1,high+(k<<20)); + return t1*w; + } else return w; +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ Index: s_logarithm.c =================================================================== --- s_logarithm.c (nonexistent) +++ s_logarithm.c (revision 1765) @@ -0,0 +1,135 @@ + +/* @(#)z_logarithm.c 1.0 98/08/13 */ +/****************************************************************** + * The following routines are coded directly from the algorithms + * and coefficients given in "Software Manual for the Elementary + * Functions" by William J. Cody, Jr. and William Waite, Prentice + * Hall, 1980. + ******************************************************************/ + +/* +FUNCTION + <>, <>, <>, <>, <>, <>---natural or base 10 logarithms + +INDEX + log +INDEX + logf +INDEX + log10 +INDEX + log10f + +ANSI_SYNOPSIS + #include + double log(double <[x]>); + float logf(float <[x]>); + double log10(double <[x]>); + float log10f(float <[x]>); + +TRAD_SYNOPSIS + #include + double log(<[x]>); + double <[x]>; + + float logf(<[x]>); + float <[x]>; + + double log10(<[x]>); + double <[x]>; + + float log10f(<[x]>); + float <[x]>; + +DESCRIPTION +Return the natural or base 10 logarithm of <[x]>, that is, its logarithm base e +(where e is the base of the natural system of logarithms, 2.71828@dots{}) or +base 10. +<> and <> are identical save for the return and argument types. +<> and <> are identical save for the return and argument types. + +RETURNS +Normally, returns the calculated value. When <[x]> is zero, the +returned value is <<-HUGE_VAL>> and <> is set to <>. +When <[x]> is negative, the returned value is <<-HUGE_VAL>> and +<> is set to <>. You can control the error behavior via +<>. + +PORTABILITY +<> is ANSI, <> is an extension. +<> is ANSI, <> is an extension. +*/ + + +/****************************************************************** + * Logarithm + * + * Input: + * x - floating point value + * ten - indicates base ten numbers + * + * Output: + * logarithm of x + * + * Description: + * This routine calculates logarithms. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +static const double a[] = { -0.64124943423745581147e+02, + 0.16383943563021534222e+02, + -0.78956112887481257267 }; +static const double b[] = { -0.76949932108494879777e+03, + 0.31203222091924532844e+03, + -0.35667977739034646171e+02 }; +static const double C1 = 22713.0 / 32768.0; +static const double C2 = 1.428606820309417232e-06; +static const double C3 = 0.43429448190325182765; + +double +_DEFUN (logarithm, (double, int), + double x _AND + int ten) +{ + int N; + double f, w, z; + + /* Check for domain error here. */ + if (x <= 0.0) + { + errno = ERANGE; + return (z_notanum.d); + } + + /* Get the exponent and mantissa where x = f * 2^N. */ + f = frexp (x, &N); + + z = f - 0.5; + + if (f > __SQRT_HALF) + z = (z - 0.5) / (f * 0.5 + 0.5); + else + { + N--; + z /= (z * 0.5 + 0.5); + } + w = z * z; + + /* Use Newton's method with 4 terms. */ + z += z * w * ((a[2] * w + a[1]) * w + a[0]) / (((w + b[2]) * w + b[1]) * w + b[0]); + + if (N != 0) + z = (N * C2 + z) + N * C1; + + if (ten) + z *= C3; + + return (z); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: sf_log.c =================================================================== --- sf_log.c (nonexistent) +++ sf_log.c (revision 1765) @@ -0,0 +1,34 @@ + +/* @(#)z_logf.c 1.0 98/08/13 */ +/****************************************************************** + * Logarithm + * + * Input: + * x - floating point value + * + * Output: + * natural logarithm of x + * + * Description: + * This routine returns the natural logarithm of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +float +_DEFUN (logf, (float), + float x) +{ + return (logarithmf (x, 0)); +} + +#ifdef _DOUBLE_IS_32BITS + +double log (double x) +{ + return (double) logf ((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */

sf_log.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: sf_tan.c =================================================================== --- sf_tan.c (nonexistent) +++ sf_tan.c (revision 1765) @@ -0,0 +1,104 @@ + +/* @(#)z_tanf.c 1.0 98/08/13 */ +/****************************************************************** + * The following routines are coded directly from the algorithms + * and coefficients given in "Software Manual for the Elementary + * Functions" by William J. Cody, Jr. and William Waite, Prentice + * Hall, 1980. + ******************************************************************/ +/****************************************************************** + * Tangent + * + * Input: + * x - floating point value + * + * Output: + * tangent of x + * + * Description: + * This routine calculates the tangent of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +static const float TWO_OVER_PI = 0.6366197723; +static const float p[] = { -0.958017723e-1 }; +static const float q[] = { -0.429135777, + 0.971685835e-2 }; + +float +_DEFUN (tanf, (float), + float x) +{ + float y, f, g, XN, xnum, xden, res; + int N; + + /* Check for special values. */ + switch (numtestf (x)) + { + case NAN: + errno = EDOM; + return (x); + case INF: + errno = EDOM; + return (z_notanum_f.f); + } + + y = fabsf (x); + + /* Check for values that are out of our range. */ + if (y > 105414357.0) + { + errno = ERANGE; + return (y); + } + + if (x < 0.0) + N = (int) (x * TWO_OVER_PI - 0.5); + else + N = (int) (x * TWO_OVER_PI + 0.5); + + XN = (float) N; + + f = x - N * __PI_OVER_TWO; + + /* Check for values that are too small. */ + if (-z_rooteps_f < f && f < z_rooteps_f) + { + xnum = f; + xden = 1.0; + } + + /* Calculate the polynomial. */ + else + { + g = f * f; + + xnum = f * (p[0] * g) + f; + xden = (q[1] * g + q[0]) * g + 1.0; + } + + /* Check for odd or even values. */ + if (N & 1) + { + xnum = -xnum; + res = xden / xnum; + } + else + { + res = xnum / xden; + } + + return (res); +} + +#ifdef _DOUBLE_IS_32BITS + +double tan (double x) +{ + return (double) tanf ((float) x); +} + +#endif /* _DOUBLE_IS_32BITS */

sf_tan.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: sf_atan.c =================================================================== --- sf_atan.c (nonexistent) +++ sf_atan.c (revision 1765) @@ -0,0 +1,45 @@ + +/* @(#)z_atanf.c 1.0 98/08/13 */ +/****************************************************************** + * Arctangent + * + * Input: + * x - floating point value + * + * Output: + * arctan of x + * + * Description: + * This routine returns the arctan of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +float +_DEFUN (atanf, (float), + float x) +{ + switch (numtestf (x)) + { + case NAN: + errno = EDOM; + return (x); + case INF: + /* this should check to see if neg NaN or pos NaN... */ + return (__PI_OVER_TWO); + case 0: + return (0.0); + default: + return (atangentf (x, 0, 0, 0)); + } +} + +#ifdef _DOUBLE_IS_32BITS +double atan (double x) +{ + return (double) atangentf ((float) x, 0, 0, 0); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */

sf_atan.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: s_asin.c =================================================================== --- s_asin.c (nonexistent) +++ s_asin.c (revision 1765) @@ -0,0 +1,29 @@ + +/* @(#)z_asin.c 1.0 98/08/13 */ +/****************************************************************** + * Arcsine + * + * Input: + * x - floating point value + * + * Output: + * arcsine of x + * + * Description: + * This routine returns the arcsine of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +double +_DEFUN (asin, (double), + double x) +{ + return (asine (x, 0)); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: sf_cos.c =================================================================== --- sf_cos.c (nonexistent) +++ sf_cos.c (revision 1765) @@ -0,0 +1,34 @@ + +/* @(#)z_cosf.c 1.0 98/08/13 */ +/****************************************************************** + * Cosine + * + * Input: + * x - floating point value + * + * Output: + * cosine of x + * + * Description: + * This routine returns the cosine of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +float +_DEFUN (cosf, (float), + float x) +{ + return (sinef (x, 1)); +} + +#ifdef _DOUBLE_IS_32BITS + +double cos (double x) +{ + return (double) sinef ((float) x, 1); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */

sf_cos.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: sf_acos.c =================================================================== --- sf_acos.c (nonexistent) +++ sf_acos.c (revision 1765) @@ -0,0 +1,33 @@ + +/* @(#)z_acosf.c 1.0 98/08/13 */ +/****************************************************************** + * Arccosine + * + * Input: + * x - floating point value + * + * Output: + * arccosine of x + * + * Description: + * This routine returns the arccosine of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +float +_DEFUN (acosf, (float), + float x) +{ + return (asinef (x, 1)); +} + +#ifdef _DOUBLE_IS_32BITS +double acos (double x) +{ + return (double) asinef ((float) x, 1); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */

sf_acos.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: ef_hypot.c =================================================================== --- ef_hypot.c (nonexistent) +++ ef_hypot.c (revision 1765) @@ -0,0 +1,82 @@ +/* ef_hypot.c -- float version of e_hypot.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ + float hypotf(float x, float y) +#else + float hypotf(x,y) + float x, y; +#endif +{ + float a=x,b=y,t1,t2,y1,y2,w; + __int32_t j,k,ha,hb; + + GET_FLOAT_WORD(ha,x); + ha &= 0x7fffffffL; + GET_FLOAT_WORD(hb,y); + hb &= 0x7fffffffL; + if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;} + SET_FLOAT_WORD(a,ha); /* a <- |a| */ + SET_FLOAT_WORD(b,hb); /* b <- |b| */ + if((ha-hb)>0xf000000L) {return a+b;} /* x/y > 2**30 */ + k=0; + if(ha > 0x58800000L) { /* a>2**50 */ + if(ha >= 0x7f800000L) { /* Inf or NaN */ + w = a+b; /* for sNaN */ + if(ha == 0x7f800000L) w = a; + if(hb == 0x7f800000L) w = b; + return w; + } + /* scale a and b by 2**-60 */ + ha -= 0x5d800000L; hb -= 0x5d800000L; k += 60; + SET_FLOAT_WORD(a,ha); + SET_FLOAT_WORD(b,hb); + } + if(hb < 0x26800000L) { /* b < 2**-50 */ + if(hb <= 0x007fffffL) { /* subnormal b or 0 */ + if(hb==0) return a; + SET_FLOAT_WORD(t1,0x3f000000L); /* t1=2^126 */ + b *= t1; + a *= t1; + k -= 126; + } else { /* scale a and b by 2^60 */ + ha += 0x5d800000; /* a *= 2^60 */ + hb += 0x5d800000; /* b *= 2^60 */ + k -= 60; + SET_FLOAT_WORD(a,ha); + SET_FLOAT_WORD(b,hb); + } + } + /* medium size a and b */ + w = a-b; + if (w>b) { + SET_FLOAT_WORD(t1,ha&0xfffff000L); + t2 = a-t1; + w = sqrtf(t1*t1-(b*(-b)-t2*(a+t1))); + } else { + a = a+a; + SET_FLOAT_WORD(y1,hb&0xfffff000L); + y2 = b - y1; + SET_FLOAT_WORD(t1,ha+0x00800000L); + t2 = a - t1; + w = sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b))); + } + if(k!=0) { + SET_FLOAT_WORD(t1,0x3f800000L+(k<<23)); + return t1*w; + } else return w; +} Index: s_cosh.c =================================================================== --- s_cosh.c (nonexistent) +++ s_cosh.c (revision 1765) @@ -0,0 +1,80 @@ + +/* @(#)z_cosh.c 1.0 98/08/13 */ + +/* + +FUNCTION + <>, <>---hyperbolic cosine + +ANSI_SYNOPSIS + #include + double cosh(double <[x]>); + float coshf(float <[x]>) + +TRAD_SYNOPSIS + #include + double cosh(<[x]>) + double <[x]>; + + float coshf(<[x]>) + float <[x]>; + +DESCRIPTION + + <> computes the hyperbolic cosine of the argument <[x]>. + <)>> is defined as + @ifinfo + . (exp(x) + exp(-x))/2 + @end ifinfo + @tex + $${(e^x + e^{-x})} \over 2$$ + @end tex + + Angles are specified in radians. + + <> is identical, save that it takes and returns <>. + +RETURNS + The computed value is returned. When the correct value would create + an overflow, <> returns the value <> with the + appropriate sign, and the global value <> is set to <>. + + You can modify error handling for these functions using the + function <>. + +PORTABILITY + <> is ANSI. + <> is an extension. + +QUICKREF + cosh ansi pure + coshf - pure +*/ + +/****************************************************************** + * Hyperbolic Cosine + * + * Input: + * x - floating point value + * + * Output: + * hyperbolic cosine of x + * + * Description: + * This routine returns the hyperbolic cosine of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +double +_DEFUN (cosh, (double), + double x) +{ + return (sineh (x, 1)); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: sf_sin.c =================================================================== --- sf_sin.c (nonexistent) +++ sf_sin.c (revision 1765) @@ -0,0 +1,34 @@ + +/* @(#)z_sinf.c 1.0 98/08/13 */ +/****************************************************************** + * Sine + * + * Input: + * x - floating point value + * + * Output: + * sine of x + * + * Description: + * This routine returns the sine of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +float +_DEFUN (sinf, (float), + float x) +{ + return (sinef (x, 0)); +} + +#ifdef _DOUBLE_IS_32BITS + +double sin (double x) +{ + return (double) sinef ((float) x, 0); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */

sf_sin.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: sf_asin.c =================================================================== --- sf_asin.c (nonexistent) +++ sf_asin.c (revision 1765) @@ -0,0 +1,34 @@ + +/* @(#)z_asinf.c 1.0 98/08/13 */ +/****************************************************************** + * Arcsine + * + * Input: + * x - floating point value + * + * Output: + * arcsine of x + * + * Description: + * This routine returns the arcsine of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +float +_DEFUN (asinf, (float), + float x) +{ + return (asinef (x, 0)); +} + +#ifdef _DOUBLE_IS_32BITS + +double asin (double x) +{ + return (double) asinef ((float) x, 0); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */

sf_asin.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: s_asinh.c =================================================================== --- s_asinh.c (nonexistent) +++ s_asinh.c (revision 1765) @@ -0,0 +1,107 @@ + +/* @(#)s_asinh.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* +FUNCTION + <>, <>---inverse hyperbolic sine + +INDEX + asinh +INDEX + asinhf + +ANSI_SYNOPSIS + #include + double asinh(double <[x]>); + float asinhf(float <[x]>); + +TRAD_SYNOPSIS + #include + double asinh(<[x]>) + double <[x]>; + + float asinhf(<[x]>) + float <[x]>; + +DESCRIPTION +<> calculates the inverse hyperbolic sine of <[x]>. +<> is defined as +@ifinfo +. sgn(<[x]>) * log(abs(<[x]>) + sqrt(1+<[x]>*<[x]>)) +@end ifinfo +@tex +$$sign(x) \times ln\Bigl(|x| + \sqrt{1+x^2}\Bigr)$$ +@end tex + +<> is identical, other than taking and returning floats. + +RETURNS +<> and <> return the calculated value. + +PORTABILITY +Neither <> nor <> are ANSI C. + +*/ + +/* asinh(x) + * Method : + * Based on + * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ] + * we have + * asinh(x) := x if 1+x*x=1, + * := sign(x)*(log(x)+ln2)) for large |x|, else + * := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else + * := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2))) + */ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ +static const double +#else +static double +#endif +one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ +ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ +huge= 1.00000000000000000000e+300; + +#ifdef __STDC__ + double asinh(double x) +#else + double asinh(x) + double x; +#endif +{ + double t,w; + __int32_t hx,ix; + GET_HIGH_WORD(hx,x); + ix = hx&0x7fffffff; + if(ix>=0x7ff00000) return x+x; /* x is inf or NaN */ + if(ix< 0x3e300000) { /* |x|<2**-28 */ + if(huge+x>one) return x; /* return x inexact except 0 */ + } + if(ix>0x41b00000) { /* |x| > 2**28 */ + w = log(fabs(x))+ln2; + } else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */ + t = fabs(x); + w = log(2.0*t+one/(sqrt(x*x+one)+t)); + } else { /* 2.0 > |x| > 2**-28 */ + t = x*x; + w =log1p(fabs(x)+t/(one+sqrt(one+t))); + } + if(hx>0) return w; else return -w; +} + +#endif /* _DOUBLE_IS_32BITS */ Index: sf_cosh.c =================================================================== --- sf_cosh.c (nonexistent) +++ sf_cosh.c (revision 1765) @@ -0,0 +1,33 @@ + +/* @(#)z_coshf.c 1.0 98/08/13 */ +/****************************************************************** + * Hyperbolic Cosine + * + * Input: + * x - floating point value + * + * Output: + * hyperbolic cosine of x + * + * Description: + * This routine returns the hyperbolic cosine of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +float +_DEFUN (coshf, (float), + float x) +{ + return (sinehf (x, 1)); +} + +#ifdef _DOUBLE_IS_32BITS +double cosh (double x) +{ + return (double) sinehf ((float) x, 1); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */

sf_cosh.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: sf_ispos.c =================================================================== --- sf_ispos.c (nonexistent) +++ sf_ispos.c (revision 1765) @@ -0,0 +1,40 @@ + +/* @(#)z_isposf.c 1.0 98/08/13 */ +/****************************************************************** + * Positive value test + * + * Input: + * x - floating point value + * + * Output: + * An integer that indicates if the number is positive. + * + * Description: + * This routine returns an integer that indicates if the number + * passed in is positive (1) or negative (0). + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +int isposf (float x) +{ + __int32_t wx; + + GET_FLOAT_WORD (wx, x); + + if (wx & 0x80000000) + return (0); + else + return (1); +} + +#ifdef _DOUBLE_IS_32BITS + +int ispos (double x) +{ + return isposf ((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */

sf_ispos.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: sf_asine.c =================================================================== --- sf_asine.c (nonexistent) +++ sf_asine.c (revision 1765) @@ -0,0 +1,105 @@ + +/* @(#)z_asinef.c 1.0 98/08/13 */ +/****************************************************************** + * The following routines are coded directly from the algorithms + * and coefficients given in "Software Manual for the Elementary + * Functions" by William J. Cody, Jr. and William Waite, Prentice + * Hall, 1980. + ******************************************************************/ +/****************************************************************** + * Arcsine + * + * Input: + * x - floating point value + * acosine - indicates acos calculation + * + * Output: + * Arcsine of x. + * + * Description: + * This routine calculates arcsine / arccosine. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +static const float p[] = { 0.933935835, -0.504400557 }; +static const float q[] = { 0.560363004e+1, -0.554846723e+1 }; +static const float a[] = { 0.0, 0.785398163 }; +static const float b[] = { 1.570796326, 0.785398163 }; + +float +_DEFUN (asinef, (float, int), + float x _AND + int acosine) +{ + int flag, i; + int branch = 0; + float g, res, R, P, Q, y; + + /* Check for special values. */ + i = numtestf (x); + if (i == NAN || i == INF) + { + errno = EDOM; + if (i == NAN) + return (x); + else + return (z_infinity_f.f); + } + + y = fabsf (x); + flag = acosine; + + if (y > 0.5) + { + i = 1 - flag; + + /* Check for range error. */ + if (y > 1.0) + { + errno = ERANGE; + return (z_notanum_f.f); + } + + g = (1 - y) / 2.0; + y = -2 * sqrt (g); + branch = 1; + } + else + { + i = flag; + if (y < z_rooteps_f) + res = y; + else + g = y * y; + } + + if (y >= z_rooteps_f || branch == 1) + { + /* Calculate the Taylor series. */ + P = (p[1] * g + p[0]) * g; + Q = (g + q[1]) * g + q[0]; + R = P / Q; + + res = y + y * R; + } + + /* Calculate asine or acose. */ + if (flag == 0) + { + res = (a[i] + res) + a[i]; + if (x < 0.0) + res = -res; + } + else + { + if (x < 0.0) + res = (b[i] + res) + b[i]; + else + res = (a[i] - res) + a[i]; + } + + return (res); +}

sf_asine.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: s_isinf.c =================================================================== --- s_isinf.c (nonexistent) +++ s_isinf.c (revision 1765) @@ -0,0 +1,37 @@ + +/* @(#)z_isinf.c 1.0 98/08/13 */ +/****************************************************************** + * isinf + * + * Input: + * x - pointer to a floating point value + * + * Output: + * An integer that indicates if the number is infinite. + * + * Description: + * This routine returns an integer that indicates if the number + * passed in is infinite (1) or is finite (0). + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +int isinf (double x) +{ + __uint32_t lx, hx; + int exp; + + EXTRACT_WORDS (hx, lx, x); + exp = (hx & 0x7ff00000) >> 20; + + if ((exp == 0x7ff) && ((hx & 0xf0000 || lx) == 0)) + return (1); + else + return (0); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: mathfp.tex =================================================================== --- mathfp.tex (nonexistent) +++ mathfp.tex (revision 1765) @@ -0,0 +1,199 @@ +@node Math +@chapter Mathematical Functions (@file{math.h}) + +This chapter groups a wide variety of mathematical functions. The +corresponding definitions and declarations are in @file{math.h}. +Two definitions from @file{math.h} are of particular interest. + +@enumerate +@item +The representation of infinity as a @code{double} is defined as +@code{HUGE_VAL}; this number is returned on overflow by many functions. + +@item +The structure @code{exception} is used when you write customized error +handlers for the mathematical functions. You can customize error +handling for most of these functions by defining your own version of +@code{matherr}; see the section on @code{matherr} for details. +@end enumerate + +@cindex system calls +@cindex support subroutines +@cindex stubs +@cindex OS stubs +Since the error handling code calls @code{fputs}, the mathematical +subroutines require stubs or minimal implementations for the same list +of OS subroutines as @code{fputs}: @code{close}, @code{fstat}, +@code{isatty}, @code{lseek}, @code{read}, @code{sbrk}, @code{write}. +@xref{syscalls,,System Calls, libc.info, The Cygnus C Support Library}, +for a discussion and for sample minimal implementations of these support +subroutines. + +Alternative declarations of the mathematical functions, which exploit +specific machine capabilities to operate faster---but generally have +less error checking and may reflect additional limitations on some +machines---are available when you include @file{fastmath.h} instead of +@file{math.h}. + +@menu +* version:: Version of library +* acos:: Arccosine +* acosh:: Inverse hyperbolic cosine +* asin:: Arcsine +* asinh:: Inverse hyperbolic sine +* atan:: Arctangent +* atan2:: Arctangent of y/x +* atanh:: Inverse hyperbolic tangent +* jN:: Bessel functions (jN, yN) +* cbrt:: Cube root +* copysign:: Sign of Y, magnitude of X +* cosh:: Hyperbolic cosine +* erf:: Error function (erf, erfc) +* exp:: Exponential +* expm1:: Exponential of x, - 1 +* fabs:: Absolute value (magnitude) +* floor:: Floor and ceiling (floor, ceil) +* fmod:: Floating-point remainder (modulo) +* frexp:: Split floating-point number +* gamma:: Logarithmic gamma function +* hypot:: Distance from origin +* ilogb:: Get exponent +* infinity:: Floating infinity +* isnan:: Check type of number +* ldexp:: Load exponent +* log:: Natural logarithms +* log10:: Base 10 logarithms +* log1p:: Log of 1 + X +* matherr:: Modifiable math error handler +* modf:: Split fractional and integer parts +* nan:: Floating Not a Number +* nextafter:: Get next representable number +* pow:: X to the power Y +* remainder:: remainder of X divided by Y +* scalbn:: scalbn +* sin:: Sine or cosine (sin, cos) +* sinh:: Hyperbolic sine +* sqrt:: Positive square root +* tan:: Tangent +* tanh:: Hyperbolic tangent +@end menu + +@page +@node version +@section Version of library + +There are four different versions of the math library routines: IEEE, +POSIX, X/Open, or SVID. The version may be selected at runtime by +setting the global variable @code{_LIB_VERSION}, defined in +@file{math.h}. It may be set to one of the following constants defined +in @file{math.h}: @code{_IEEE_}, @code{_POSIX_}, @code{_XOPEN_}, or +@code{_SVID_}. The @code{_LIB_VERSION} variable is not specific to any +thread, and changing it will affect all threads. + +The versions of the library differ only in how errors are handled. + +In IEEE mode, the @code{matherr} function is never called, no warning +messages are printed, and @code{errno} is never set. + +In POSIX mode, @code{errno} is set correctly, but the @code{matherr} +function is never called and no warning messages are printed. + +In X/Open mode, @code{errno} is set correctly, and @code{matherr} is +called, but warning message are not printed. + +In SVID mode, functions which overflow return 3.40282346638528860e+38, +the maximum single precision floating point value, rather than infinity. +Also, @code{errno} is set correctly, @code{matherr} is called, and, if +@code{matherr} returns 0, warning messages are printed for some errors. +For example, by default @samp{log(-1.0)} writes this message on standard +error output: + +@example +log: DOMAIN error +@end example + +The library is set to X/Open mode by default. + +@page +@include mathfp/sacos.def + +@page +@include mathfp/eacosh.def + +@page +@include mathfp/sasine.def + +@page +@include mathfp/sasinh.def + +@page +@include mathfp/satan.def + +@page +@include mathfp/satan2.def + +@page +@include mathfp/eatanh.def + +@page +@include mathfp/wjn.def + +@page +@include mathfp/scosh.def + +@page +@include mathfp/serf.def + +@page +@include mathfp/sexp.def + +@page +@include mathfp/sfabs.def + +@page +@include mathfp/sfloor.def + +@page +@include mathfp/sfmod.def + +@page +@include mathfp/sfrexp.def + +@page +@include mathfp/erlgamma.def + +@page +@include mathfp/ehypot.def + +@page +@include mathfp/sisnan.def + +@page +@include mathfp/sldexp.def + +@page +@include mathfp/slogarithm.def + +@page +@include mathfp/slog10.def + +@page +@include mathfp/spow.def + +@page +@include mathfp/eremainder.def + +@page +@include mathfp/ssqrt.def + +@page +@include mathfp/ssine.def + +@page +@include mathfp/ssineh.def + +@page +@include mathfp/stan.def + +@page +@include mathfp/stanh.def Index: sf_pow.c =================================================================== --- sf_pow.c (nonexistent) +++ sf_pow.c (revision 1765) @@ -0,0 +1,107 @@ + +/* @(#)z_powf.c 1.0 98/08/13 */ +#include +#include "fdlibm.h" +#include "zmath.h" + +float powf (float x, float y) +{ + float d, t, r = 1.0; + int n, k, sign = 0; + __int32_t px; + + GET_FLOAT_WORD (px, x); + + k = modff (y, &d); + if (k == 0.0) + { + if (modff (ldexpf (y, -1), &t)) + sign = 0; + else + sign = 1; + } + + if (x == 0.0 && y <= 0.0) + errno = EDOM; + + else if ((t = y * log (fabsf (x))) >= BIGX) + { + errno = ERANGE; + if (px & 0x80000000) + { + if (!k) + { + errno = EDOM; + x = 0.0; + } + else if (sign) + x = -z_infinity_f.f; + else + x = z_infinity_f.f; + } + + else + x = z_infinity_f.f; + } + + else if (t < SMALLX) + { + errno = ERANGE; + x = 0.0; + } + + else + { + if ( k && fabsf (d) <= 32767 ) + { + n = (int) d; + + if (sign = (n < 0)) + n = -n; + + while ( n > 0 ) + { + if ((unsigned int) n % 2) + r *= x; + x *= x; + n = (unsigned int) n / 2; + } + + if (sign) + r = 1.0 / r; + + return r; + } + + else + { + if ( px & 0x80000000 ) + { + if ( !k ) + { + errno = EDOM; + return 0.0; + } + } + + x = exp (t); + + if ( sign ) + { + px ^= 0x80000000; + SET_FLOAT_WORD (x, px); + } + } + } + + return x; +} + +#ifdef _DOUBLE_IS_32BITS + +double pow (double x, double y) +{ + return (double) powf ((float) x, (float) y); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ Index: s_erf.c =================================================================== --- s_erf.c (nonexistent) +++ s_erf.c (revision 1765) @@ -0,0 +1,373 @@ + +/* @(#)s_erf.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* +FUNCTION + <>, <>, <>, <>---error function +INDEX + erf +INDEX + erff +INDEX + erfc +INDEX + erfcf + +ANSI_SYNOPSIS + #include + double erf(double <[x]>); + float erff(float <[x]>); + double erfc(double <[x]>); + float erfcf(float <[x]>); +TRAD_SYNOPSIS + #include + + double erf(<[x]>) + double <[x]>; + + float erff(<[x]>) + float <[x]>; + + double erfc(<[x]>) + double <[x]>; + + float erfcf(<[x]>) + float <[x]>; + +DESCRIPTION + <> calculates an approximation to the ``error function'', + which estimates the probability that an observation will fall within + <[x]> standard deviations of the mean (assuming a normal + distribution). + @tex + The error function is defined as + $${2\over\sqrt\pi}\times\int_0^x e^{-t^2}dt$$ + @end tex + + <> calculates the complementary probability; that is, + <)>> is <<1 - erf(<[x]>)>>. <> is computed directly, + so that you can use it to avoid the loss of precision that would + result from subtracting large probabilities (on large <[x]>) from 1. + + <> and <> differ from <> and <> only in the + argument and result types. + +RETURNS + For positive arguments, <> and all its variants return a + probability---a number between 0 and 1. + +PORTABILITY + None of the variants of <> are ANSI C. +*/ + +/* double erf(double x) + * double erfc(double x) + * x + * 2 |\ + * erf(x) = --------- | exp(-t*t)dt + * sqrt(pi) \| + * 0 + * + * erfc(x) = 1-erf(x) + * Note that + * erf(-x) = -erf(x) + * erfc(-x) = 2 - erfc(x) + * + * Method: + * 1. For |x| in [0, 0.84375] + * erf(x) = x + x*R(x^2) + * erfc(x) = 1 - erf(x) if x in [-.84375,0.25] + * = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375] + * where R = P/Q where P is an odd poly of degree 8 and + * Q is an odd poly of degree 10. + * -57.90 + * | R - (erf(x)-x)/x | <= 2 + * + * + * Remark. The formula is derived by noting + * erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....) + * and that + * 2/sqrt(pi) = 1.128379167095512573896158903121545171688 + * is close to one. The interval is chosen because the fix + * point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is + * near 0.6174), and by some experiment, 0.84375 is chosen to + * guarantee the error is less than one ulp for erf. + * + * 2. For |x| in [0.84375,1.25], let s = |x| - 1, and + * c = 0.84506291151 rounded to single (24 bits) + * erf(x) = sign(x) * (c + P1(s)/Q1(s)) + * erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0 + * 1+(c+P1(s)/Q1(s)) if x < 0 + * |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06 + * Remark: here we use the taylor series expansion at x=1. + * erf(1+s) = erf(1) + s*Poly(s) + * = 0.845.. + P1(s)/Q1(s) + * That is, we use rational approximation to approximate + * erf(1+s) - (c = (single)0.84506291151) + * Note that |P1/Q1|< 0.078 for x in [0.84375,1.25] + * where + * P1(s) = degree 6 poly in s + * Q1(s) = degree 6 poly in s + * + * 3. For x in [1.25,1/0.35(~2.857143)], + * erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1) + * erf(x) = 1 - erfc(x) + * where + * R1(z) = degree 7 poly in z, (z=1/x^2) + * S1(z) = degree 8 poly in z + * + * 4. For x in [1/0.35,28] + * erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0 + * = 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6 x >= 28 + * erf(x) = sign(x) *(1 - tiny) (raise inexact) + * erfc(x) = tiny*tiny (raise underflow) if x > 0 + * = 2 - tiny if x<0 + * + * 7. Special case: + * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1, + * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2, + * erfc/erf(NaN) is NaN + */ + + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ +static const double +#else +static double +#endif +tiny = 1e-300, +half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ +one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ +two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */ + /* c = (float)0.84506291151 */ +erx = 8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */ +/* + * Coefficients for approximation to erf on [0,0.84375] + */ +efx = 1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */ +efx8= 1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */ +pp0 = 1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */ +pp1 = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */ +pp2 = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */ +pp3 = -5.77027029648944159157e-03, /* 0xBF77A291, 0x236668E4 */ +pp4 = -2.37630166566501626084e-05, /* 0xBEF8EAD6, 0x120016AC */ +qq1 = 3.97917223959155352819e-01, /* 0x3FD97779, 0xCDDADC09 */ +qq2 = 6.50222499887672944485e-02, /* 0x3FB0A54C, 0x5536CEBA */ +qq3 = 5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */ +qq4 = 1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */ +qq5 = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */ +/* + * Coefficients for approximation to erf in [0.84375,1.25] + */ +pa0 = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */ +pa1 = 4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */ +pa2 = -3.72207876035701323847e-01, /* 0xBFD7D240, 0xFBB8C3F1 */ +pa3 = 3.18346619901161753674e-01, /* 0x3FD45FCA, 0x805120E4 */ +pa4 = -1.10894694282396677476e-01, /* 0xBFBC6398, 0x3D3E28EC */ +pa5 = 3.54783043256182359371e-02, /* 0x3FA22A36, 0x599795EB */ +pa6 = -2.16637559486879084300e-03, /* 0xBF61BF38, 0x0A96073F */ +qa1 = 1.06420880400844228286e-01, /* 0x3FBB3E66, 0x18EEE323 */ +qa2 = 5.40397917702171048937e-01, /* 0x3FE14AF0, 0x92EB6F33 */ +qa3 = 7.18286544141962662868e-02, /* 0x3FB2635C, 0xD99FE9A7 */ +qa4 = 1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */ +qa5 = 1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */ +qa6 = 1.19844998467991074170e-02, /* 0x3F888B54, 0x5735151D */ +/* + * Coefficients for approximation to erfc in [1.25,1/0.35] + */ +ra0 = -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */ +ra1 = -6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */ +ra2 = -1.05586262253232909814e+01, /* 0xC0251E04, 0x41B0E726 */ +ra3 = -6.23753324503260060396e+01, /* 0xC04F300A, 0xE4CBA38D */ +ra4 = -1.62396669462573470355e+02, /* 0xC0644CB1, 0x84282266 */ +ra5 = -1.84605092906711035994e+02, /* 0xC067135C, 0xEBCCABB2 */ +ra6 = -8.12874355063065934246e+01, /* 0xC0545265, 0x57E4D2F2 */ +ra7 = -9.81432934416914548592e+00, /* 0xC023A0EF, 0xC69AC25C */ +sa1 = 1.96512716674392571292e+01, /* 0x4033A6B9, 0xBD707687 */ +sa2 = 1.37657754143519042600e+02, /* 0x4061350C, 0x526AE721 */ +sa3 = 4.34565877475229228821e+02, /* 0x407B290D, 0xD58A1A71 */ +sa4 = 6.45387271733267880336e+02, /* 0x40842B19, 0x21EC2868 */ +sa5 = 4.29008140027567833386e+02, /* 0x407AD021, 0x57700314 */ +sa6 = 1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */ +sa7 = 6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */ +sa8 = -6.04244152148580987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */ +/* + * Coefficients for approximation to erfc in [1/.35,28] + */ +rb0 = -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */ +rb1 = -7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */ +rb2 = -1.77579549177547519889e+01, /* 0xC031C209, 0x555F995A */ +rb3 = -1.60636384855821916062e+02, /* 0xC064145D, 0x43C5ED98 */ +rb4 = -6.37566443368389627722e+02, /* 0xC083EC88, 0x1375F228 */ +rb5 = -1.02509513161107724954e+03, /* 0xC0900461, 0x6A2E5992 */ +rb6 = -4.83519191608651397019e+02, /* 0xC07E384E, 0x9BDC383F */ +sb1 = 3.03380607434824582924e+01, /* 0x403E568B, 0x261D5190 */ +sb2 = 3.25792512996573918826e+02, /* 0x40745CAE, 0x221B9F0A */ +sb3 = 1.53672958608443695994e+03, /* 0x409802EB, 0x189D5118 */ +sb4 = 3.19985821950859553908e+03, /* 0x40A8FFB7, 0x688C246A */ +sb5 = 2.55305040643316442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */ +sb6 = 4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */ +sb7 = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */ + +#ifdef __STDC__ + double erf(double x) +#else + double erf(x) + double x; +#endif +{ + __int32_t hx,ix,i; + double R,S,P,Q,s,y,z,r; + GET_HIGH_WORD(hx,x); + ix = hx&0x7fffffff; + if(ix>=0x7ff00000) { /* erf(nan)=nan */ + i = ((__uint32_t)hx>>31)<<1; + return (double)(1-i)+one/x; /* erf(+-inf)=+-1 */ + } + + if(ix < 0x3feb0000) { /* |x|<0.84375 */ + if(ix < 0x3e300000) { /* |x|<2**-28 */ + if (ix < 0x00800000) + return 0.125*(8.0*x+efx8*x); /*avoid underflow */ + return x + efx*x; + } + z = x*x; + r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4))); + s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); + y = r/s; + return x + x*y; + } + if(ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */ + s = fabs(x)-one; + P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); + Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6))))); + if(hx>=0) return erx + P/Q; else return -erx - P/Q; + } + if (ix >= 0x40180000) { /* inf>|x|>=6 */ + if(hx>=0) return one-tiny; else return tiny-one; + } + x = fabs(x); + s = one/(x*x); + if(ix< 0x4006DB6E) { /* |x| < 1/0.35 */ + R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( + ra5+s*(ra6+s*ra7)))))); + S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*( + sa5+s*(sa6+s*(sa7+s*sa8))))))); + } else { /* |x| >= 1/0.35 */ + R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( + rb5+s*rb6))))); + S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( + sb5+s*(sb6+s*sb7)))))); + } + z = x; + SET_LOW_WORD(z,0); + r = exp(-z*z-0.5625)*exp((z-x)*(z+x)+R/S); + if(hx>=0) return one-r/x; else return r/x-one; +} + +#ifdef __STDC__ + double erfc(double x) +#else + double erfc(x) + double x; +#endif +{ + __int32_t hx,ix; + double R,S,P,Q,s,y,z,r; + GET_HIGH_WORD(hx,x); + ix = hx&0x7fffffff; + if(ix>=0x7ff00000) { /* erfc(nan)=nan */ + /* erfc(+-inf)=0,2 */ + return (double)(((__uint32_t)hx>>31)<<1)+one/x; + } + + if(ix < 0x3feb0000) { /* |x|<0.84375 */ + if(ix < 0x3c700000) /* |x|<2**-56 */ + return one-x; + z = x*x; + r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4))); + s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); + y = r/s; + if(hx < 0x3fd00000) { /* x<1/4 */ + return one-(x+x*y); + } else { + r = x*y; + r += (x-half); + return half - r ; + } + } + if(ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */ + s = fabs(x)-one; + P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); + Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6))))); + if(hx>=0) { + z = one-erx; return z - P/Q; + } else { + z = erx+P/Q; return one+z; + } + } + if (ix < 0x403c0000) { /* |x|<28 */ + x = fabs(x); + s = one/(x*x); + if(ix< 0x4006DB6D) { /* |x| < 1/.35 ~ 2.857143*/ + R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( + ra5+s*(ra6+s*ra7)))))); + S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*( + sa5+s*(sa6+s*(sa7+s*sa8))))))); + } else { /* |x| >= 1/.35 ~ 2.857143 */ + if(hx<0&&ix>=0x40180000) return two-tiny;/* x < -6 */ + R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( + rb5+s*rb6))))); + S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( + sb5+s*(sb6+s*sb7)))))); + } + z = x; + SET_LOW_WORD(z,0); + r = exp(-z*z-0.5625)* + exp((z-x)*(z+x)+R/S); + if(hx>0) return r/x; else return two-r/x; + } else { + if(hx>0) return tiny*tiny; else return two-tiny; + } +} + +#endif /* _DOUBLE_IS_32BITS */ Index: e_atanh.c =================================================================== --- e_atanh.c (nonexistent) +++ e_atanh.c (revision 1765) @@ -0,0 +1,139 @@ + +/* @(#)e_atanh.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + */ + +/* +FUNCTION + <>, <>---inverse hyperbolic tangent + +INDEX + atanh +INDEX + atanhf + +ANSI_SYNOPSIS + #include + double atanh(double <[x]>); + float atanhf(float <[x]>); + +TRAD_SYNOPSIS + #include + double atanh(<[x]>) + double <[x]>; + + float atanhf(<[x]>) + float <[x]>; + +DESCRIPTION + <> calculates the inverse hyperbolic tangent of <[x]>. + + <> is identical, other than taking and returning + <> values. + +RETURNS + <> and <> return the calculated value. + + If + @ifinfo + |<[x]>| + @end ifinfo + @tex + $|x|$ + @end tex + is greater than 1, the global <> is set to <> and + the result is a NaN. A <> is reported. + + If + @ifinfo + |<[x]>| + @end ifinfo + @tex + $|x|$ + @end tex + is 1, the global <> is set to <>; and the result is + infinity with the same sign as <>. A <> is reported. + + You can modify the error handling for these routines using + <>. + +PORTABILITY + Neither <> nor <> are ANSI C. + +QUICKREF + atanh - pure + atanhf - pure + + +*/ + +/* atanh(x) + * Method : + * 1.Reduced x to positive by atanh(-x) = -atanh(x) + * 2.For x>=0.5 + * 1 2x x + * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) + * 2 1 - x 1 - x + * + * For x<0.5 + * atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) + * + * Special cases: + * atanh(x) is NaN if |x| > 1 with signal; + * atanh(NaN) is that NaN with no signal; + * atanh(+-1) is +-INF with signal. + * + */ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ +static const double one = 1.0, huge = 1e300; +#else +static double one = 1.0, huge = 1e300; +#endif + +#ifdef __STDC__ +static const double zero = 0.0; +#else +static double zero = 0.0; +#endif + +#ifdef __STDC__ + double atanh(double x) +#else + double atanh(x) + double x; +#endif +{ + double t; + __int32_t hx,ix; + __uint32_t lx; + EXTRACT_WORDS(hx,lx,x); + ix = hx&0x7fffffff; + if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */ + return (x-x)/(x-x); + if(ix==0x3ff00000) + return x/zero; + if(ix<0x3e300000&&(huge+x)>zero) return x; /* x<2**-28 */ + SET_HIGH_WORD(x,ix); + if(ix<0x3fe00000) { /* x < 0.5 */ + t = x+x; + t = 0.5*log1p(t+t*x/(one-x)); + } else + t = 0.5*log1p((x+x)/(one-x)); + if(hx>=0) return t; else return -t; +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ Index: ef_atanh.c =================================================================== --- ef_atanh.c (nonexistent) +++ ef_atanh.c (revision 1765) @@ -0,0 +1,54 @@ +/* ef_atanh.c -- float version of e_atanh.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + */ + +#include "fdlibm.h" + +#ifdef __STDC__ +static const float one = 1.0, huge = 1e30; +#else +static float one = 1.0, huge = 1e30; +#endif + +#ifdef __STDC__ +static const float zero = 0.0; +#else +static float zero = 0.0; +#endif + +#ifdef __STDC__ + float atanhf(float x) +#else + float atanhf(x) + float x; +#endif +{ + float t; + __int32_t hx,ix; + GET_FLOAT_WORD(hx,x); + ix = hx&0x7fffffff; + if (ix>0x3f800000) /* |x|>1 */ + return (x-x)/(x-x); + if(ix==0x3f800000) + return x/zero; + if(ix<0x31800000&&(huge+x)>zero) return x; /* x<2**-28 */ + SET_FLOAT_WORD(x,ix); + if(ix<0x3f000000) { /* x < 0.5 */ + t = x+x; + t = (float)0.5*log1pf(t+t*x/(one-x)); + } else + t = (float)0.5*log1pf((x+x)/(one-x)); + if(hx>=0) return t; else return -t; +} Index: er_gamma.c =================================================================== --- er_gamma.c (nonexistent) +++ er_gamma.c (revision 1765) @@ -0,0 +1,32 @@ + +/* @(#)er_gamma.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + */ + +/* gamma_r(x, signgamp) + * Reentrant version of the logarithm of the Gamma function + * with user provide pointer for the sign of Gamma(x). + * + * Method: See lgamma_r + */ + +#include "fdlibm.h" + +#ifdef __STDC__ + double gamma_r(double x, int *signgamp) +#else + double gamma_r(x,signgamp) + double x; int *signgamp; +#endif +{ + return lgamma_r(x,signgamp); +} Index: sf_floor.c =================================================================== --- sf_floor.c (nonexistent) +++ sf_floor.c (revision 1765) @@ -0,0 +1,43 @@ + +/* @(#)z_floorf.c 1.0 98/08/13 */ +/***************************************************************** + * floor + * + * Input: + * x - floating point value + * + * Output: + * Smallest integer less than x. + * + * Description: + * This routine returns the smallest integer less than x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +float +_DEFUN (floorf, (float), + float x) +{ + float f, y; + + if (x > -1.0 && x < 1.0) + return (x >= 0 ? 0 : -1.0); + + y = modff (x, &f); + + if (y == 0.0) + return (x); + + return (x >= 0 ? f : f - 1.0); +} + +#ifdef _DOUBLE_IS_32BITS +double floor (double x) +{ + return (double) floorf ((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */

sf_floor.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: e_remainder.c =================================================================== --- e_remainder.c (nonexistent) +++ e_remainder.c (revision 1765) @@ -0,0 +1,113 @@ + +/* @(#)e_remainder.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* +FUNCTION +<>, <>---round and remainder +INDEX + remainder +INDEX + remainderf + +ANSI_SYNOPSIS + #include + double remainder(double <[x]>, double <[y]>); + float remainderf(float <[x]>, float <[y]>); + +TRAD_SYNOPSIS + #include + double remainder(<[x]>,<[y]>) + double <[x]>, <[y]>; + float remainderf(<[x]>,<[y]>) + float <[x]>, <[y]>; + +DESCRIPTION +<> and <> find the remainder of +<[x]>/<[y]>; this value is in the range -<[y]>/2 .. +<[y]>/2. + +RETURNS +<> returns the integer result as a double. + +PORTABILITY +<> is a System V release 4. +<> is an extension. + +*/ + +/* remainder(x,p) + * Return : + * returns x REM p = x - [x/p]*p as if in infinite + * precise arithmetic, where [x/p] is the (infinite bit) + * integer nearest x/p (in half way case choose the even one). + * Method : + * Based on fmod() return x-[x/p]chopped*p exactlp. + */ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ +static const double zero = 0.0; +#else +static double zero = 0.0; +#endif + + +#ifdef __STDC__ + double remainder(double x, double p) +#else + double remainder(x,p) + double x,p; +#endif +{ + __int32_t hx,hp; + __uint32_t sx,lx,lp; + double p_half; + + EXTRACT_WORDS(hx,lx,x); + EXTRACT_WORDS(hp,lp,p); + sx = hx&0x80000000; + hp &= 0x7fffffff; + hx &= 0x7fffffff; + + /* purge off exception values */ + if((hp|lp)==0) return (x*p)/(x*p); /* p = 0 */ + if((hx>=0x7ff00000)|| /* x not finite */ + ((hp>=0x7ff00000)&& /* p is NaN */ + (((hp-0x7ff00000)|lp)!=0))) + return (x*p)/(x*p); + + + if (hp<=0x7fdfffff) x = fmod(x,p+p); /* now x < 2p */ + if (((hx-hp)|(lx-lp))==0) return zero*x; + x = fabs(x); + p = fabs(p); + if (hp<0x00200000) { + if(x+x>p) { + x-=p; + if(x+x>=p) x -= p; + } + } else { + p_half = 0.5*p; + if(x>p_half) { + x-=p; + if(x>=p_half) x -= p; + } + } + GET_HIGH_WORD(hx,x); + SET_HIGH_WORD(x,hx^sx); + return x; +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ Index: w_drem.c =================================================================== --- w_drem.c (nonexistent) +++ w_drem.c (revision 1765) @@ -0,0 +1,15 @@ +/* + * drem() wrapper for remainder(). + * + * Written by J.T. Conklin, + * Placed into the Public Domain, 1994. + */ + +#include "fdlibm.h" + +double +drem(x, y) + double x, y; +{ + return remainder(x, y); +} Index: e_j1.c =================================================================== --- e_j1.c (nonexistent) +++ e_j1.c (revision 1765) @@ -0,0 +1,486 @@ + +/* @(#)e_j1.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* j1(x), y1(x) + * Bessel function of the first and second kinds of order zero. + * Method -- j1(x): + * 1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ... + * 2. Reduce x to |x| since j1(x)=-j1(-x), and + * for x in (0,2) + * j1(x) = x/2 + x*z*R0/S0, where z = x*x; + * (precision: |j1/x - 1/2 - R0/S0 |<2**-61.51 ) + * for x in (2,inf) + * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1)) + * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) + * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) + * as follow: + * cos(x1) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) + * = 1/sqrt(2) * (sin(x) - cos(x)) + * sin(x1) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) + * = -1/sqrt(2) * (sin(x) + cos(x)) + * (To avoid cancellation, use + * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) + * to compute the worse one.) + * + * 3 Special cases + * j1(nan)= nan + * j1(0) = 0 + * j1(inf) = 0 + * + * Method -- y1(x): + * 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN + * 2. For x<2. + * Since + * y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...) + * therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function. + * We use the following function to approximate y1, + * y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2 + * where for x in [0,2] (abs err less than 2**-65.89) + * U(z) = U0[0] + U0[1]*z + ... + U0[4]*z^4 + * V(z) = 1 + v0[0]*z + ... + v0[4]*z^5 + * Note: For tiny x, 1/x dominate y1 and hence + * y1(tiny) = -2/pi/tiny, (choose tiny<2**-54) + * 3. For x>=2. + * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1)) + * where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1) + * by method mentioned above. + */ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +#ifdef __STDC__ +static double pone(double), qone(double); +#else +static double pone(), qone(); +#endif + +#ifdef __STDC__ +static const double +#else +static double +#endif +huge = 1e300, +one = 1.0, +invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */ +tpi = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ + /* R0/S0 on [0,2] */ +r00 = -6.25000000000000000000e-02, /* 0xBFB00000, 0x00000000 */ +r01 = 1.40705666955189706048e-03, /* 0x3F570D9F, 0x98472C61 */ +r02 = -1.59955631084035597520e-05, /* 0xBEF0C5C6, 0xBA169668 */ +r03 = 4.96727999609584448412e-08, /* 0x3E6AAAFA, 0x46CA0BD9 */ +s01 = 1.91537599538363460805e-02, /* 0x3F939D0B, 0x12637E53 */ +s02 = 1.85946785588630915560e-04, /* 0x3F285F56, 0xB9CDF664 */ +s03 = 1.17718464042623683263e-06, /* 0x3EB3BFF8, 0x333F8498 */ +s04 = 5.04636257076217042715e-09, /* 0x3E35AC88, 0xC97DFF2C */ +s05 = 1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */ + +#ifdef __STDC__ +static const double zero = 0.0; +#else +static double zero = 0.0; +#endif + +#ifdef __STDC__ + double j1(double x) +#else + double j1(x) + double x; +#endif +{ + double z, s,c,ss,cc,r,u,v,y; + __int32_t hx,ix; + + GET_HIGH_WORD(hx,x); + ix = hx&0x7fffffff; + if(ix>=0x7ff00000) return one/x; + y = fabs(x); + if(ix >= 0x40000000) { /* |x| >= 2.0 */ + s = sin(y); + c = cos(y); + ss = -s-c; + cc = s-c; + if(ix<0x7fe00000) { /* make sure y+y not overflow */ + z = cos(y+y); + if ((s*c)>zero) cc = z/ss; + else ss = z/cc; + } + /* + * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x) + * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x) + */ + if(ix>0x48000000) z = (invsqrtpi*cc)/sqrt(y); + else { + u = pone(y); v = qone(y); + z = invsqrtpi*(u*cc-v*ss)/sqrt(y); + } + if(hx<0) return -z; + else return z; + } + if(ix<0x3e400000) { /* |x|<2**-27 */ + if(huge+x>one) return 0.5*x;/* inexact if x!=0 necessary */ + } + z = x*x; + r = z*(r00+z*(r01+z*(r02+z*r03))); + s = one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05)))); + r *= x; + return(x*0.5+r/s); +} + +#ifdef __STDC__ +static const double U0[5] = { +#else +static double U0[5] = { +#endif + -1.96057090646238940668e-01, /* 0xBFC91866, 0x143CBC8A */ + 5.04438716639811282616e-02, /* 0x3FA9D3C7, 0x76292CD1 */ + -1.91256895875763547298e-03, /* 0xBF5F55E5, 0x4844F50F */ + 2.35252600561610495928e-05, /* 0x3EF8AB03, 0x8FA6B88E */ + -9.19099158039878874504e-08, /* 0xBE78AC00, 0x569105B8 */ +}; +#ifdef __STDC__ +static const double V0[5] = { +#else +static double V0[5] = { +#endif + 1.99167318236649903973e-02, /* 0x3F94650D, 0x3F4DA9F0 */ + 2.02552581025135171496e-04, /* 0x3F2A8C89, 0x6C257764 */ + 1.35608801097516229404e-06, /* 0x3EB6C05A, 0x894E8CA6 */ + 6.22741452364621501295e-09, /* 0x3E3ABF1D, 0x5BA69A86 */ + 1.66559246207992079114e-11, /* 0x3DB25039, 0xDACA772A */ +}; + +#ifdef __STDC__ + double y1(double x) +#else + double y1(x) + double x; +#endif +{ + double z, s,c,ss,cc,u,v; + __int32_t hx,ix,lx; + + EXTRACT_WORDS(hx,lx,x); + ix = 0x7fffffff&hx; + /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ + if(ix>=0x7ff00000) return one/(x+x*x); + if((ix|lx)==0) return -one/zero; + if(hx<0) return zero/zero; + if(ix >= 0x40000000) { /* |x| >= 2.0 */ + s = sin(x); + c = cos(x); + ss = -s-c; + cc = s-c; + if(ix<0x7fe00000) { /* make sure x+x not overflow */ + z = cos(x+x); + if ((s*c)>zero) cc = z/ss; + else ss = z/cc; + } + /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0)) + * where x0 = x-3pi/4 + * Better formula: + * cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4) + * = 1/sqrt(2) * (sin(x) - cos(x)) + * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) + * = -1/sqrt(2) * (cos(x) + sin(x)) + * To avoid cancellation, use + * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) + * to compute the worse one. + */ + if(ix>0x48000000) z = (invsqrtpi*ss)/sqrt(x); + else { + u = pone(x); v = qone(x); + z = invsqrtpi*(u*ss+v*cc)/sqrt(x); + } + return z; + } + if(ix<=0x3c900000) { /* x < 2**-54 */ + return(-tpi/x); + } + z = x*x; + u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4]))); + v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4])))); + return(x*(u/v) + tpi*(j1(x)*log(x)-one/x)); +} + +/* For x >= 8, the asymptotic expansions of pone is + * 1 + 15/128 s^2 - 4725/2^15 s^4 - ..., where s = 1/x. + * We approximate pone by + * pone(x) = 1 + (R/S) + * where R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10 + * S = 1 + ps0*s^2 + ... + ps4*s^10 + * and + * | pone(x)-1-R/S | <= 2 ** ( -60.06) + */ + +#ifdef __STDC__ +static const double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ +#else +static double pr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ +#endif + 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ + 1.17187499999988647970e-01, /* 0x3FBDFFFF, 0xFFFFFCCE */ + 1.32394806593073575129e+01, /* 0x402A7A9D, 0x357F7FCE */ + 4.12051854307378562225e+02, /* 0x4079C0D4, 0x652EA590 */ + 3.87474538913960532227e+03, /* 0x40AE457D, 0xA3A532CC */ + 7.91447954031891731574e+03, /* 0x40BEEA7A, 0xC32782DD */ +}; +#ifdef __STDC__ +static const double ps8[5] = { +#else +static double ps8[5] = { +#endif + 1.14207370375678408436e+02, /* 0x405C8D45, 0x8E656CAC */ + 3.65093083420853463394e+03, /* 0x40AC85DC, 0x964D274F */ + 3.69562060269033463555e+04, /* 0x40E20B86, 0x97C5BB7F */ + 9.76027935934950801311e+04, /* 0x40F7D42C, 0xB28F17BB */ + 3.08042720627888811578e+04, /* 0x40DE1511, 0x697A0B2D */ +}; + +#ifdef __STDC__ +static const double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ +#else +static double pr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ +#endif + 1.31990519556243522749e-11, /* 0x3DAD0667, 0xDAE1CA7D */ + 1.17187493190614097638e-01, /* 0x3FBDFFFF, 0xE2C10043 */ + 6.80275127868432871736e+00, /* 0x401B3604, 0x6E6315E3 */ + 1.08308182990189109773e+02, /* 0x405B13B9, 0x452602ED */ + 5.17636139533199752805e+02, /* 0x40802D16, 0xD052D649 */ + 5.28715201363337541807e+02, /* 0x408085B8, 0xBB7E0CB7 */ +}; +#ifdef __STDC__ +static const double ps5[5] = { +#else +static double ps5[5] = { +#endif + 5.92805987221131331921e+01, /* 0x404DA3EA, 0xA8AF633D */ + 9.91401418733614377743e+02, /* 0x408EFB36, 0x1B066701 */ + 5.35326695291487976647e+03, /* 0x40B4E944, 0x5706B6FB */ + 7.84469031749551231769e+03, /* 0x40BEA4B0, 0xB8A5BB15 */ + 1.50404688810361062679e+03, /* 0x40978030, 0x036F5E51 */ +}; + +#ifdef __STDC__ +static const double pr3[6] = { +#else +static double pr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ +#endif + 3.02503916137373618024e-09, /* 0x3E29FC21, 0xA7AD9EDD */ + 1.17186865567253592491e-01, /* 0x3FBDFFF5, 0x5B21D17B */ + 3.93297750033315640650e+00, /* 0x400F76BC, 0xE85EAD8A */ + 3.51194035591636932736e+01, /* 0x40418F48, 0x9DA6D129 */ + 9.10550110750781271918e+01, /* 0x4056C385, 0x4D2C1837 */ + 4.85590685197364919645e+01, /* 0x4048478F, 0x8EA83EE5 */ +}; +#ifdef __STDC__ +static const double ps3[5] = { +#else +static double ps3[5] = { +#endif + 3.47913095001251519989e+01, /* 0x40416549, 0xA134069C */ + 3.36762458747825746741e+02, /* 0x40750C33, 0x07F1A75F */ + 1.04687139975775130551e+03, /* 0x40905B7C, 0x5037D523 */ + 8.90811346398256432622e+02, /* 0x408BD67D, 0xA32E31E9 */ + 1.03787932439639277504e+02, /* 0x4059F26D, 0x7C2EED53 */ +}; + +#ifdef __STDC__ +static const double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ +#else +static double pr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ +#endif + 1.07710830106873743082e-07, /* 0x3E7CE9D4, 0xF65544F4 */ + 1.17176219462683348094e-01, /* 0x3FBDFF42, 0xBE760D83 */ + 2.36851496667608785174e+00, /* 0x4002F2B7, 0xF98FAEC0 */ + 1.22426109148261232917e+01, /* 0x40287C37, 0x7F71A964 */ + 1.76939711271687727390e+01, /* 0x4031B1A8, 0x177F8EE2 */ + 5.07352312588818499250e+00, /* 0x40144B49, 0xA574C1FE */ +}; +#ifdef __STDC__ +static const double ps2[5] = { +#else +static double ps2[5] = { +#endif + 2.14364859363821409488e+01, /* 0x40356FBD, 0x8AD5ECDC */ + 1.25290227168402751090e+02, /* 0x405F5293, 0x14F92CD5 */ + 2.32276469057162813669e+02, /* 0x406D08D8, 0xD5A2DBD9 */ + 1.17679373287147100768e+02, /* 0x405D6B7A, 0xDA1884A9 */ + 8.36463893371618283368e+00, /* 0x4020BAB1, 0xF44E5192 */ +}; + +#ifdef __STDC__ + static double pone(double x) +#else + static double pone(x) + double x; +#endif +{ +#ifdef __STDC__ + const double *p,*q; +#else + double *p,*q; +#endif + double z,r,s; + __int32_t ix; + GET_HIGH_WORD(ix,x); + ix &= 0x7fffffff; + if(ix>=0x40200000) {p = pr8; q= ps8;} + else if(ix>=0x40122E8B){p = pr5; q= ps5;} + else if(ix>=0x4006DB6D){p = pr3; q= ps3;} + else {p = pr2; q= ps2;} + z = one/(x*x); + r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); + s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); + return one+ r/s; +} + + +/* For x >= 8, the asymptotic expansions of qone is + * 3/8 s - 105/1024 s^3 - ..., where s = 1/x. + * We approximate qone by + * qone(x) = s*(0.375 + (R/S)) + * where R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10 + * S = 1 + qs1*s^2 + ... + qs6*s^12 + * and + * | qone(x)/s -0.375-R/S | <= 2 ** ( -61.13) + */ + +#ifdef __STDC__ +static const double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ +#else +static double qr8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ +#endif + 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ + -1.02539062499992714161e-01, /* 0xBFBA3FFF, 0xFFFFFDF3 */ + -1.62717534544589987888e+01, /* 0xC0304591, 0xA26779F7 */ + -7.59601722513950107896e+02, /* 0xC087BCD0, 0x53E4B576 */ + -1.18498066702429587167e+04, /* 0xC0C724E7, 0x40F87415 */ + -4.84385124285750353010e+04, /* 0xC0E7A6D0, 0x65D09C6A */ +}; +#ifdef __STDC__ +static const double qs8[6] = { +#else +static double qs8[6] = { +#endif + 1.61395369700722909556e+02, /* 0x40642CA6, 0xDE5BCDE5 */ + 7.82538599923348465381e+03, /* 0x40BE9162, 0xD0D88419 */ + 1.33875336287249578163e+05, /* 0x4100579A, 0xB0B75E98 */ + 7.19657723683240939863e+05, /* 0x4125F653, 0x72869C19 */ + 6.66601232617776375264e+05, /* 0x412457D2, 0x7719AD5C */ + -2.94490264303834643215e+05, /* 0xC111F969, 0x0EA5AA18 */ +}; + +#ifdef __STDC__ +static const double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ +#else +static double qr5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ +#endif + -2.08979931141764104297e-11, /* 0xBDB6FA43, 0x1AA1A098 */ + -1.02539050241375426231e-01, /* 0xBFBA3FFF, 0xCB597FEF */ + -8.05644828123936029840e+00, /* 0xC0201CE6, 0xCA03AD4B */ + -1.83669607474888380239e+02, /* 0xC066F56D, 0x6CA7B9B0 */ + -1.37319376065508163265e+03, /* 0xC09574C6, 0x6931734F */ + -2.61244440453215656817e+03, /* 0xC0A468E3, 0x88FDA79D */ +}; +#ifdef __STDC__ +static const double qs5[6] = { +#else +static double qs5[6] = { +#endif + 8.12765501384335777857e+01, /* 0x405451B2, 0xFF5A11B2 */ + 1.99179873460485964642e+03, /* 0x409F1F31, 0xE77BF839 */ + 1.74684851924908907677e+04, /* 0x40D10F1F, 0x0D64CE29 */ + 4.98514270910352279316e+04, /* 0x40E8576D, 0xAABAD197 */ + 2.79480751638918118260e+04, /* 0x40DB4B04, 0xCF7C364B */ + -4.71918354795128470869e+03, /* 0xC0B26F2E, 0xFCFFA004 */ +}; + +#ifdef __STDC__ +static const double qr3[6] = { +#else +static double qr3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ +#endif + -5.07831226461766561369e-09, /* 0xBE35CFA9, 0xD38FC84F */ + -1.02537829820837089745e-01, /* 0xBFBA3FEB, 0x51AEED54 */ + -4.61011581139473403113e+00, /* 0xC01270C2, 0x3302D9FF */ + -5.78472216562783643212e+01, /* 0xC04CEC71, 0xC25D16DA */ + -2.28244540737631695038e+02, /* 0xC06C87D3, 0x4718D55F */ + -2.19210128478909325622e+02, /* 0xC06B66B9, 0x5F5C1BF6 */ +}; +#ifdef __STDC__ +static const double qs3[6] = { +#else +static double qs3[6] = { +#endif + 4.76651550323729509273e+01, /* 0x4047D523, 0xCCD367E4 */ + 6.73865112676699709482e+02, /* 0x40850EEB, 0xC031EE3E */ + 3.38015286679526343505e+03, /* 0x40AA684E, 0x448E7C9A */ + 5.54772909720722782367e+03, /* 0x40B5ABBA, 0xA61D54A6 */ + 1.90311919338810798763e+03, /* 0x409DBC7A, 0x0DD4DF4B */ + -1.35201191444307340817e+02, /* 0xC060E670, 0x290A311F */ +}; + +#ifdef __STDC__ +static const double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ +#else +static double qr2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ +#endif + -1.78381727510958865572e-07, /* 0xBE87F126, 0x44C626D2 */ + -1.02517042607985553460e-01, /* 0xBFBA3E8E, 0x9148B010 */ + -2.75220568278187460720e+00, /* 0xC0060484, 0x69BB4EDA */ + -1.96636162643703720221e+01, /* 0xC033A9E2, 0xC168907F */ + -4.23253133372830490089e+01, /* 0xC04529A3, 0xDE104AAA */ + -2.13719211703704061733e+01, /* 0xC0355F36, 0x39CF6E52 */ +}; +#ifdef __STDC__ +static const double qs2[6] = { +#else +static double qs2[6] = { +#endif + 2.95333629060523854548e+01, /* 0x403D888A, 0x78AE64FF */ + 2.52981549982190529136e+02, /* 0x406F9F68, 0xDB821CBA */ + 7.57502834868645436472e+02, /* 0x4087AC05, 0xCE49A0F7 */ + 7.39393205320467245656e+02, /* 0x40871B25, 0x48D4C029 */ + 1.55949003336666123687e+02, /* 0x40637E5E, 0x3C3ED8D4 */ + -4.95949898822628210127e+00, /* 0xC013D686, 0xE71BE86B */ +}; + +#ifdef __STDC__ + static double qone(double x) +#else + static double qone(x) + double x; +#endif +{ +#ifdef __STDC__ + const double *p,*q; +#else + double *p,*q; +#endif + double s,r,z; + __int32_t ix; + GET_HIGH_WORD(ix,x); + ix &= 0x7fffffff; + if(ix>=0x40200000) {p = qr8; q= qs8;} + else if(ix>=0x40122E8B){p = qr5; q= qs5;} + else if(ix>=0x4006DB6D){p = qr3; q= qs3;} + else {p = qr2; q= qs2;} + z = one/(x*x); + r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); + s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); + return (.375 + r/s)/x; +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ Index: ef_remainder.c =================================================================== --- ef_remainder.c (nonexistent) +++ ef_remainder.c (revision 1765) @@ -0,0 +1,68 @@ +/* ef_remainder.c -- float version of e_remainder.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ +static const float zero = 0.0; +#else +static float zero = 0.0; +#endif + + +#ifdef __STDC__ + float remainderf(float x, float p) +#else + float remainderf(x,p) + float x,p; +#endif +{ + __int32_t hx,hp; + __uint32_t sx; + float p_half; + + GET_FLOAT_WORD(hx,x); + GET_FLOAT_WORD(hp,p); + sx = hx&0x80000000; + hp &= 0x7fffffff; + hx &= 0x7fffffff; + + /* purge off exception values */ + if(hp==0) return (x*p)/(x*p); /* p = 0 */ + if((hx>=0x7f800000)|| /* x not finite */ + ((hp>0x7f800000))) /* p is NaN */ + return (x*p)/(x*p); + + + if (hp<=0x7effffff) x = fmodf(x,p+p); /* now x < 2p */ + if ((hx-hp)==0) return zero*x; + x = fabsf(x); + p = fabsf(p); + if (hp<0x01000000) { + if(x+x>p) { + x-=p; + if(x+x>=p) x -= p; + } + } else { + p_half = (float)0.5*p; + if(x>p_half) { + x-=p; + if(x>=p_half) x -= p; + } + } + GET_FLOAT_WORD(hx,x); + SET_FLOAT_WORD(x,hx^sx); + return x; +} Index: s_sinf.c =================================================================== --- s_sinf.c (nonexistent) +++ s_sinf.c (revision 1765) @@ -0,0 +1,34 @@ + +/* @(#)z_sinf.c 1.0 98/08/13 */ +/****************************************************************** + * Sine + * + * Input: + * x - floating point value + * + * Output: + * sine of x + * + * Description: + * This routine returns the sine of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +float +_DEFUN (sinf, (float), + float x) +{ + return (sinef (x, 0)); +} + +#ifdef _DOUBLE_IS_32BITS + +double sin (double x) +{ + return (double) sinf ((float) x); +} + +#endif /* _DOUBLE_IS_32BITS */

s_sinf.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: s_sinh.c =================================================================== --- s_sinh.c (nonexistent) +++ s_sinh.c (revision 1765) @@ -0,0 +1,29 @@ + +/* @(#)z_sinh.c 1.0 98/08/13 */ +/****************************************************************** + * Hyperbolic Sine + * + * Input: + * x - floating point value + * + * Output: + * hyperbolic sine of x + * + * Description: + * This routine returns the hyperbolic sine of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +double +_DEFUN (sinh, (double), + double x) +{ + return (sineh (x, 0)); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: ef_j0.c =================================================================== --- ef_j0.c (nonexistent) +++ ef_j0.c (revision 1765) @@ -0,0 +1,439 @@ +/* ef_j0.c -- float version of e_j0.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ +static float pzerof(float), qzerof(float); +#else +static float pzerof(), qzerof(); +#endif + +#ifdef __STDC__ +static const float +#else +static float +#endif +huge = 1e30, +one = 1.0, +invsqrtpi= 5.6418961287e-01, /* 0x3f106ebb */ +tpi = 6.3661974669e-01, /* 0x3f22f983 */ + /* R0/S0 on [0, 2.00] */ +R02 = 1.5625000000e-02, /* 0x3c800000 */ +R03 = -1.8997929874e-04, /* 0xb947352e */ +R04 = 1.8295404516e-06, /* 0x35f58e88 */ +R05 = -4.6183270541e-09, /* 0xb19eaf3c */ +S01 = 1.5619102865e-02, /* 0x3c7fe744 */ +S02 = 1.1692678527e-04, /* 0x38f53697 */ +S03 = 5.1354652442e-07, /* 0x3509daa6 */ +S04 = 1.1661400734e-09; /* 0x30a045e8 */ + +#ifdef __STDC__ +static const float zero = 0.0; +#else +static float zero = 0.0; +#endif + +#ifdef __STDC__ + float j0f(float x) +#else + float j0f(x) + float x; +#endif +{ + float z, s,c,ss,cc,r,u,v; + __int32_t hx,ix; + + GET_FLOAT_WORD(hx,x); + ix = hx&0x7fffffff; + if(ix>=0x7f800000) return one/(x*x); + x = fabsf(x); + if(ix >= 0x40000000) { /* |x| >= 2.0 */ + s = sinf(x); + c = cosf(x); + ss = s-c; + cc = s+c; + if(ix<0x7f000000) { /* make sure x+x not overflow */ + z = -cosf(x+x); + if ((s*c)0x80000000) z = (invsqrtpi*cc)/sqrtf(x); + else { + u = pzerof(x); v = qzerof(x); + z = invsqrtpi*(u*cc-v*ss)/sqrtf(x); + } + return z; + } + if(ix<0x39000000) { /* |x| < 2**-13 */ + if(huge+x>one) { /* raise inexact if x != 0 */ + if(ix<0x32000000) return one; /* |x|<2**-27 */ + else return one - (float)0.25*x*x; + } + } + z = x*x; + r = z*(R02+z*(R03+z*(R04+z*R05))); + s = one+z*(S01+z*(S02+z*(S03+z*S04))); + if(ix < 0x3F800000) { /* |x| < 1.00 */ + return one + z*((float)-0.25+(r/s)); + } else { + u = (float)0.5*x; + return((one+u)*(one-u)+z*(r/s)); + } +} + +#ifdef __STDC__ +static const float +#else +static float +#endif +u00 = -7.3804296553e-02, /* 0xbd9726b5 */ +u01 = 1.7666645348e-01, /* 0x3e34e80d */ +u02 = -1.3818567619e-02, /* 0xbc626746 */ +u03 = 3.4745343146e-04, /* 0x39b62a69 */ +u04 = -3.8140706238e-06, /* 0xb67ff53c */ +u05 = 1.9559013964e-08, /* 0x32a802ba */ +u06 = -3.9820518410e-11, /* 0xae2f21eb */ +v01 = 1.2730483897e-02, /* 0x3c509385 */ +v02 = 7.6006865129e-05, /* 0x389f65e0 */ +v03 = 2.5915085189e-07, /* 0x348b216c */ +v04 = 4.4111031494e-10; /* 0x2ff280c2 */ + +#ifdef __STDC__ + float y0f(float x) +#else + float y0f(x) + float x; +#endif +{ + float z, s,c,ss,cc,u,v; + __int32_t hx,ix; + + GET_FLOAT_WORD(hx,x); + ix = 0x7fffffff&hx; + /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */ + if(ix>=0x7f800000) return one/(x+x*x); + if(ix==0) return -one/zero; + if(hx<0) return zero/zero; + if(ix >= 0x40000000) { /* |x| >= 2.0 */ + /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0)) + * where x0 = x-pi/4 + * Better formula: + * cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4) + * = 1/sqrt(2) * (sin(x) + cos(x)) + * sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4) + * = 1/sqrt(2) * (sin(x) - cos(x)) + * To avoid cancellation, use + * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) + * to compute the worse one. + */ + s = sinf(x); + c = cosf(x); + ss = s-c; + cc = s+c; + /* + * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x) + * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x) + */ + if(ix<0x7f000000) { /* make sure x+x not overflow */ + z = -cosf(x+x); + if ((s*c)0x80000000) z = (invsqrtpi*ss)/sqrtf(x); + else { + u = pzerof(x); v = qzerof(x); + z = invsqrtpi*(u*ss+v*cc)/sqrtf(x); + } + return z; + } + if(ix<=0x32000000) { /* x < 2**-27 */ + return(u00 + tpi*logf(x)); + } + z = x*x; + u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06))))); + v = one+z*(v01+z*(v02+z*(v03+z*v04))); + return(u/v + tpi*(j0f(x)*logf(x))); +} + +/* The asymptotic expansions of pzero is + * 1 - 9/128 s^2 + 11025/98304 s^4 - ..., where s = 1/x. + * For x >= 2, We approximate pzero by + * pzero(x) = 1 + (R/S) + * where R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10 + * S = 1 + pS0*s^2 + ... + pS4*s^10 + * and + * | pzero(x)-1-R/S | <= 2 ** ( -60.26) + */ +#ifdef __STDC__ +static const float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ +#else +static float pR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ +#endif + 0.0000000000e+00, /* 0x00000000 */ + -7.0312500000e-02, /* 0xbd900000 */ + -8.0816707611e+00, /* 0xc1014e86 */ + -2.5706311035e+02, /* 0xc3808814 */ + -2.4852163086e+03, /* 0xc51b5376 */ + -5.2530439453e+03, /* 0xc5a4285a */ +}; +#ifdef __STDC__ +static const float pS8[5] = { +#else +static float pS8[5] = { +#endif + 1.1653436279e+02, /* 0x42e91198 */ + 3.8337448730e+03, /* 0x456f9beb */ + 4.0597855469e+04, /* 0x471e95db */ + 1.1675296875e+05, /* 0x47e4087c */ + 4.7627726562e+04, /* 0x473a0bba */ +}; +#ifdef __STDC__ +static const float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ +#else +static float pR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ +#endif + -1.1412546255e-11, /* 0xad48c58a */ + -7.0312492549e-02, /* 0xbd8fffff */ + -4.1596107483e+00, /* 0xc0851b88 */ + -6.7674766541e+01, /* 0xc287597b */ + -3.3123129272e+02, /* 0xc3a59d9b */ + -3.4643338013e+02, /* 0xc3ad3779 */ +}; +#ifdef __STDC__ +static const float pS5[5] = { +#else +static float pS5[5] = { +#endif + 6.0753936768e+01, /* 0x42730408 */ + 1.0512523193e+03, /* 0x44836813 */ + 5.9789707031e+03, /* 0x45bad7c4 */ + 9.6254453125e+03, /* 0x461665c8 */ + 2.4060581055e+03, /* 0x451660ee */ +}; + +#ifdef __STDC__ +static const float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ +#else +static float pR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ +#endif + -2.5470459075e-09, /* 0xb12f081b */ + -7.0311963558e-02, /* 0xbd8fffb8 */ + -2.4090321064e+00, /* 0xc01a2d95 */ + -2.1965976715e+01, /* 0xc1afba52 */ + -5.8079170227e+01, /* 0xc2685112 */ + -3.1447946548e+01, /* 0xc1fb9565 */ +}; +#ifdef __STDC__ +static const float pS3[5] = { +#else +static float pS3[5] = { +#endif + 3.5856033325e+01, /* 0x420f6c94 */ + 3.6151397705e+02, /* 0x43b4c1ca */ + 1.1936077881e+03, /* 0x44953373 */ + 1.1279968262e+03, /* 0x448cffe6 */ + 1.7358093262e+02, /* 0x432d94b8 */ +}; + +#ifdef __STDC__ +static const float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ +#else +static float pR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ +#endif + -8.8753431271e-08, /* 0xb3be98b7 */ + -7.0303097367e-02, /* 0xbd8ffb12 */ + -1.4507384300e+00, /* 0xbfb9b1cc */ + -7.6356959343e+00, /* 0xc0f4579f */ + -1.1193166733e+01, /* 0xc1331736 */ + -3.2336456776e+00, /* 0xc04ef40d */ +}; +#ifdef __STDC__ +static const float pS2[5] = { +#else +static float pS2[5] = { +#endif + 2.2220300674e+01, /* 0x41b1c32d */ + 1.3620678711e+02, /* 0x430834f0 */ + 2.7047027588e+02, /* 0x43873c32 */ + 1.5387539673e+02, /* 0x4319e01a */ + 1.4657617569e+01, /* 0x416a859a */ +}; + +#ifdef __STDC__ + static float pzerof(float x) +#else + static float pzerof(x) + float x; +#endif +{ +#ifdef __STDC__ + const float *p,*q; +#else + float *p,*q; +#endif + float z,r,s; + __int32_t ix; + GET_FLOAT_WORD(ix,x); + ix &= 0x7fffffff; + if(ix>=0x41000000) {p = pR8; q= pS8;} + else if(ix>=0x40f71c58){p = pR5; q= pS5;} + else if(ix>=0x4036db68){p = pR3; q= pS3;} + else {p = pR2; q= pS2;} + z = one/(x*x); + r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); + s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); + return one+ r/s; +} + + +/* For x >= 8, the asymptotic expansions of qzero is + * -1/8 s + 75/1024 s^3 - ..., where s = 1/x. + * We approximate qzero by + * qzero(x) = s*(-1.25 + (R/S)) + * where R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10 + * S = 1 + qS0*s^2 + ... + qS5*s^12 + * and + * | qzero(x)/s +1.25-R/S | <= 2 ** ( -61.22) + */ +#ifdef __STDC__ +static const float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ +#else +static float qR8[6] = { /* for x in [inf, 8]=1/[0,0.125] */ +#endif + 0.0000000000e+00, /* 0x00000000 */ + 7.3242187500e-02, /* 0x3d960000 */ + 1.1768206596e+01, /* 0x413c4a93 */ + 5.5767340088e+02, /* 0x440b6b19 */ + 8.8591972656e+03, /* 0x460a6cca */ + 3.7014625000e+04, /* 0x471096a0 */ +}; +#ifdef __STDC__ +static const float qS8[6] = { +#else +static float qS8[6] = { +#endif + 1.6377603149e+02, /* 0x4323c6aa */ + 8.0983447266e+03, /* 0x45fd12c2 */ + 1.4253829688e+05, /* 0x480b3293 */ + 8.0330925000e+05, /* 0x49441ed4 */ + 8.4050156250e+05, /* 0x494d3359 */ + -3.4389928125e+05, /* 0xc8a7eb69 */ +}; + +#ifdef __STDC__ +static const float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ +#else +static float qR5[6] = { /* for x in [8,4.5454]=1/[0.125,0.22001] */ +#endif + 1.8408595828e-11, /* 0x2da1ec79 */ + 7.3242180049e-02, /* 0x3d95ffff */ + 5.8356351852e+00, /* 0x40babd86 */ + 1.3511157227e+02, /* 0x43071c90 */ + 1.0272437744e+03, /* 0x448067cd */ + 1.9899779053e+03, /* 0x44f8bf4b */ +}; +#ifdef __STDC__ +static const float qS5[6] = { +#else +static float qS5[6] = { +#endif + 8.2776611328e+01, /* 0x42a58da0 */ + 2.0778142090e+03, /* 0x4501dd07 */ + 1.8847289062e+04, /* 0x46933e94 */ + 5.6751113281e+04, /* 0x475daf1d */ + 3.5976753906e+04, /* 0x470c88c1 */ + -5.3543427734e+03, /* 0xc5a752be */ +}; + +#ifdef __STDC__ +static const float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ +#else +static float qR3[6] = {/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */ +#endif + 4.3774099900e-09, /* 0x3196681b */ + 7.3241114616e-02, /* 0x3d95ff70 */ + 3.3442313671e+00, /* 0x405607e3 */ + 4.2621845245e+01, /* 0x422a7cc5 */ + 1.7080809021e+02, /* 0x432acedf */ + 1.6673394775e+02, /* 0x4326bbe4 */ +}; +#ifdef __STDC__ +static const float qS3[6] = { +#else +static float qS3[6] = { +#endif + 4.8758872986e+01, /* 0x42430916 */ + 7.0968920898e+02, /* 0x44316c1c */ + 3.7041481934e+03, /* 0x4567825f */ + 6.4604252930e+03, /* 0x45c9e367 */ + 2.5163337402e+03, /* 0x451d4557 */ + -1.4924745178e+02, /* 0xc3153f59 */ +}; + +#ifdef __STDC__ +static const float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ +#else +static float qR2[6] = {/* for x in [2.8570,2]=1/[0.3499,0.5] */ +#endif + 1.5044444979e-07, /* 0x342189db */ + 7.3223426938e-02, /* 0x3d95f62a */ + 1.9981917143e+00, /* 0x3fffc4bf */ + 1.4495602608e+01, /* 0x4167edfd */ + 3.1666231155e+01, /* 0x41fd5471 */ + 1.6252708435e+01, /* 0x4182058c */ +}; +#ifdef __STDC__ +static const float qS2[6] = { +#else +static float qS2[6] = { +#endif + 3.0365585327e+01, /* 0x41f2ecb8 */ + 2.6934811401e+02, /* 0x4386ac8f */ + 8.4478375244e+02, /* 0x44533229 */ + 8.8293585205e+02, /* 0x445cbbe5 */ + 2.1266638184e+02, /* 0x4354aa98 */ + -5.3109550476e+00, /* 0xc0a9f358 */ +}; + +#ifdef __STDC__ + static float qzerof(float x) +#else + static float qzerof(x) + float x; +#endif +{ +#ifdef __STDC__ + const float *p,*q; +#else + float *p,*q; +#endif + float s,r,z; + __int32_t ix; + GET_FLOAT_WORD(ix,x); + ix &= 0x7fffffff; + if(ix>=0x41000000) {p = qR8; q= qS8;} + else if(ix>=0x40f71c58){p = qR5; q= qS5;} + else if(ix>=0x4036db68){p = qR3; q= qS3;} + else {p = qR2; q= qS2;} + z = one/(x*x); + r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5])))); + s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5]))))); + return (-(float).125 + r/s)/x; +} Index: wf_drem.c =================================================================== --- wf_drem.c (nonexistent) +++ wf_drem.c (revision 1765) @@ -0,0 +1,19 @@ +/* + * dremf() wrapper for remainderf(). + * + * Written by J.T. Conklin, + * Placed into the Public Domain, 1994. + */ + +#include "fdlibm.h" + +float +#ifdef __STDC__ +dremf(float x, float y) +#else +dremf(x, y) + float x, y; +#endif +{ + return remainderf(x, y); +} Index: sf_exp.c =================================================================== --- sf_exp.c (nonexistent) +++ sf_exp.c (revision 1765) @@ -0,0 +1,92 @@ + +/* @(#)z_expf.c 1.0 98/08/13 */ +/****************************************************************** + * The following routines are coded directly from the algorithms + * and coefficients given in "Software Manual for the Elementary + * Functions" by William J. Cody, Jr. and William Waite, Prentice + * Hall, 1980. + ******************************************************************/ +/****************************************************************** + * Exponential Function + * + * Input: + * x - floating point value + * + * Output: + * e raised to x. + * + * Description: + * This routine returns e raised to the xth power. + * + *****************************************************************/ + +#include +#include "fdlibm.h" +#include "zmath.h" + +static const float INV_LN2 = 1.442695040; +static const float LN2 = 0.693147180; +static const float p[] = { 0.249999999950, 0.00416028863 }; +static const float q[] = { 0.5, 0.04998717878 }; + +float +_DEFUN (expf, (float), + float x) +{ + int N; + float g, z, R, P, Q; + + switch (numtestf (x)) + { + case NAN: + errno = EDOM; + return (x); + case INF: + errno = ERANGE; + if (isposf (x)) + return (z_infinity_f.f); + else + return (0.0); + case 0: + return (1.0); + } + + /* Check for out of bounds. */ + if (x > BIGX || x < SMALLX) + { + errno = ERANGE; + return (x); + } + + /* Check for a value too small to calculate. */ + if (-z_rooteps_f < x && x < z_rooteps_f) + { + return (1.0); + } + + /* Calculate the exponent. */ + if (x < 0.0) + N = (int) (x * INV_LN2 - 0.5); + else + N = (int) (x * INV_LN2 + 0.5); + + /* Construct the mantissa. */ + g = x - N * LN2; + z = g * g; + P = g * (p[1] * z + p[0]); + Q = q[1] * z + q[0]; + R = 0.5 + P / (Q - P); + + /* Return the floating point value. */ + N++; + return (ldexpf (R, N)); +} + +#ifdef _DOUBLE_IS_32BITS + +double exp (double x) +{ + return (double) expf ((float) x); +} + +#endif /* _DOUBLE_IS_32BITS */

sf_exp.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: sf_numtest.c =================================================================== --- sf_numtest.c (nonexistent) +++ sf_numtest.c (revision 1765) @@ -0,0 +1,63 @@ + +/* @(#)z_numtestf.c 1.0 98/08/13 */ +/****************************************************************** + * Numtest + * + * Input: + * x - pointer to a floating point value + * + * Output: + * An integer that indicates what kind of number was passed in: + * NUM = 3 - a finite value + * NAN = 2 - not a number + * INF = 1 - an infinite value + * 0 - zero + * + * Description: + * This routine returns an integer that indicates the character- + * istics of the number that was passed in. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +int +_DEFUN (numtestf, (float), + float x) +{ + __int32_t wx; + int exp; + + GET_FLOAT_WORD (wx, x); + + exp = (wx & 0x7f800000) >> 23; + + /* Check for a zero input. */ + if (x == 0.0) + { + return (0); + } + + /* Check for not a number or infinity. */ + if (exp == 0x7f8) + { + if(wx & 0x7fffff) + return (NAN); + else + return (INF); + } + + /* Otherwise it's a finite value. */ + else + return (NUM); +} + +#ifdef _DOUBLE_IS_32BITS + +int numtest (double x) +{ + return numtestf ((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */

sf_numtest.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: e_scalb.c =================================================================== --- e_scalb.c (nonexistent) +++ e_scalb.c (revision 1765) @@ -0,0 +1,55 @@ + +/* @(#)e_scalb.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + * __ieee754_scalb(x, fn) is provide for + * passing various standard test suite. One + * should use scalbn() instead. + */ + +#include "fdlibm.h" + +#ifndef _DOUBLE_IS_32BITS + +#ifdef _SCALB_INT +#ifdef __STDC__ + double scalb(double x, int fn) +#else + double scalb(x,fn) + double x; int fn; +#endif +#else +#ifdef __STDC__ + double scalb(double x, double fn) +#else + double scalb(x,fn) + double x, fn; +#endif +#endif +{ +#ifdef _SCALB_INT + return scalbn(x,fn); +#else + if (isnan(x)||isnan(fn)) return x*fn; + if (!finite(fn)) { + if(fn>0.0) return x*fn; + else return x/(-fn); + } + if (rint(fn)!=fn) return (fn-fn)/(fn-fn); + if ( fn > 65000.0) return scalbn(x, 65000); + if (-fn > 65000.0) return scalbn(x,-65000); + return scalbn(x,(int)fn); +#endif +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ Index: sf_sinh.c =================================================================== --- sf_sinh.c (nonexistent) +++ sf_sinh.c (revision 1765) @@ -0,0 +1,34 @@ + +/* @(#)z_sinhf.c 1.0 98/08/13 */ +/****************************************************************** + * Hyperbolic Sine + * + * Input: + * x - floating point value + * + * Output: + * hyperbolic sine of x + * + * Description: + * This routine returns the hyperbolic sine of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +float +_DEFUN (sinhf, (float), + float x) +{ + return (sinehf (x, 0)); +} + +#ifdef _DOUBLE_IS_32BITS + +double sinh (double x) +{ + return (double) sinhf ((float) x); +} + +#endif /* _DOUBLE_IS_32BITS */

sf_sinh.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: s_isnan.c =================================================================== --- s_isnan.c (nonexistent) +++ s_isnan.c (revision 1765) @@ -0,0 +1,125 @@ + +/* @(#)z_isnan.c 1.0 98/08/13 */ + +/* +FUNCTION + <>,<>,<>,<>,<>,<>---test +for exceptional numbers + +INDEX + isnan +INDEX + isinf +INDEX + finite + +INDEX + isnanf +INDEX + isinff +INDEX + finitef + +ANSI_SYNOPSIS + #include + int isnan(double <[arg]>); + int isinf(double <[arg]>); + int finite(double <[arg]>); + int isnanf(float <[arg]>); + int isinff(float <[arg]>); + int finitef(float <[arg]>); + +TRAD_SYNOPSIS + #include + int isnan(<[arg]>) + double <[arg]>; + int isinf(<[arg]>) + double <[arg]>; + int finite(<[arg]>); + double <[arg]>; + int isnanf(<[arg]>); + float <[arg]>; + int isinff(<[arg]>); + float <[arg]>; + int finitef(<[arg]>); + float <[arg]>; + + +DESCRIPTION + These functions provide information on the floating point + argument supplied. + + There are five major number formats - + o+ + o zero + a number which contains all zero bits. + o subnormal + Is used to represent number with a zero exponent, but a non zero fract +ion. + o normal + A number with an exponent, and a fraction + o infinity + A number with an all 1's exponent and a zero fraction. + o NAN + A number with an all 1's exponent and a non zero fraction. + + o- + + <> returns 1 if the argument is a nan. <> + returns 1 if the argument is infinity. <> returns 1 if the + argument is zero, subnormal or normal. + + The <>, <> and <> perform the same + operations as their <>, <> and <> + counterparts, but on single precision floating point numbers. + +QUICKREF + isnan - pure +QUICKREF + isinf - pure +QUICKREF + finite - pure +QUICKREF + isnan - pure +QUICKREF + isinf - pure +QUICKREF + finite - pure +*/ + + +/****************************************************************** + * isnan + * + * Input: + * x - pointer to a floating point value + * + * Output: + * An integer that indicates if the number is NaN. + * + * Description: + * This routine returns an integer that indicates if the number + * passed in is NaN (1) or is finite (0). + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +int isnan (double x) +{ + __uint32_t lx, hx; + int exp; + + EXTRACT_WORDS (hx, lx, x); + exp = (hx & 0x7ff00000) >> 20; + + if ((exp == 0x7ff) && (hx & 0xf0000 || lx)) + return (1); + else + return (0); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: sf_atan2.c =================================================================== --- sf_atan2.c (nonexistent) +++ sf_atan2.c (revision 1765) @@ -0,0 +1,34 @@ + +/* @(#)z_atan2f.c 1.0 98/08/13 */ +/****************************************************************** + * Arctangent2 + * + * Input: + * v, u - floating point values + * + * Output: + * arctan2 of v / u + * + * Description: + * This routine returns the arctan2 of v / u. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +float +_DEFUN (atan2f, (float, float), + float v _AND + float u) +{ + return (atangentf (0.0, v, u, 1)); +} + +#ifdef _DOUBLE_IS_32BITS +double atan2 (double v, double u) +{ + return (double) atangentf (0.0, (float) v, (float) u, 1); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */

sf_atan2.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: s_fabs.c =================================================================== --- s_fabs.c (nonexistent) +++ s_fabs.c (revision 1765) @@ -0,0 +1,80 @@ + +/* @(#)z_fabs.c 1.0 98/08/13 */ + +/* +FUNCTION + <>, <>---absolute value (magnitude) +INDEX + fabs +INDEX + fabsf + +ANSI_SYNOPSIS + #include + double fabs(double <[x]>); + float fabsf(float <[x]>); + +TRAD_SYNOPSIS + #include + double fabs(<[x]>) + double <[x]>; + + float fabsf(<[x]>) + float <[x]>; + +DESCRIPTION +<> and <> calculate +@tex +$|x|$, +@end tex +the absolute value (magnitude) of the argument <[x]>, by direct +manipulation of the bit representation of <[x]>. + +RETURNS +The calculated value is returned. + +PORTABILITY +<> is ANSI. +<> is an extension. + +*/ + +/****************************************************************** + * Floating-Point Absolute Value + * + * Input: + * x - floating-point number + * + * Output: + * absolute value of x + * + * Description: + * fabs computes the absolute value of a floating point number. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +double +_DEFUN (fabs, (double), + double x) +{ + switch (numtest (x)) + { + case NAN: + errno = EDOM; + return (x); + case INF: + errno = ERANGE; + return (x); + case 0: + return (0.0); + default: + return (x < 0.0 ? -x : x); + } +} + +#endif /* _DOUBLE_IS_32BITS */ Index: s_ldexp.c =================================================================== --- s_ldexp.c (nonexistent) +++ s_ldexp.c (revision 1765) @@ -0,0 +1,125 @@ + +/* @(#)z_ldexp.c 1.0 98/08/13 */ + +/* +FUNCTION + <>, <>---load exponent + +INDEX + ldexp +INDEX + ldexpf + +ANSI_SYNOPSIS + #include + double ldexp(double <[val]>, int <[exp]>); + float ldexpf(float <[val]>, int <[exp]>); + +TRAD_SYNOPSIS + #include + + double ldexp(<[val]>, <[exp]>) + double <[val]>; + int <[exp]>; + + float ldexpf(<[val]>, <[exp]>) + float <[val]>; + int <[exp]>; + +DESCRIPTION +<> calculates the value +@ifinfo +<[val]> times 2 to the power <[exp]>. +@end ifinfo +@tex +$val\times 2^{exp}$. +@end tex +<> is identical, save that it takes and returns <> +rather than <> values. + +RETURNS +<> returns the calculated value. + +Underflow and overflow both set <> to <>. +On underflow, <> and <> return 0.0. +On overflow, <> returns plus or minus <>. + +PORTABILITY +<> is ANSI, <> is an extension. + +*/ + +/****************************************************************** + * ldexp + * + * Input: + * d - a floating point value + * e - an exponent value + * + * Output: + * A floating point value f such that f = d * 2 ^ e. + * + * Description: + * This function creates a floating point number f such that + * f = d * 2 ^ e. + * + *****************************************************************/ + +#include +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +#define DOUBLE_EXP_OFFS 1023 + +double +_DEFUN (ldexp, (double, int), + double d _AND + int e) +{ + int exp; + __uint32_t hd; + + GET_HIGH_WORD (hd, d); + + /* Check for special values and then scale d by e. */ + switch (numtest (d)) + { + case NAN: + errno = EDOM; + break; + + case INF: + errno = ERANGE; + break; + + case 0: + break; + + default: + exp = (hd & 0x7ff00000) >> 20; + exp += e; + + if (exp > DBL_MAX_EXP + DOUBLE_EXP_OFFS) + { + errno = ERANGE; + d = z_infinity.d; + } + else if (exp < DBL_MIN_EXP + DOUBLE_EXP_OFFS) + { + errno = ERANGE; + d = -z_infinity.d; + } + else + { + hd &= 0x800fffff; + hd |= exp << 20; + SET_HIGH_WORD (d, hd); + } + } + + return (d); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: ef_scalb.c =================================================================== --- ef_scalb.c (nonexistent) +++ ef_scalb.c (revision 1765) @@ -0,0 +1,53 @@ +/* ef_scalb.c -- float version of e_scalb.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" +#include + +#ifdef _SCALB_INT +#ifdef __STDC__ + float scalbf(float x, int fn) +#else + float scalbf(x,fn) + float x; int fn; +#endif +#else +#ifdef __STDC__ + float scalbf(float x, float fn) +#else + float scalbf(x,fn) + float x, fn; +#endif +#endif +{ +#ifdef _SCALB_INT + return scalbnf(x,fn); +#else + if (isnanf(x)||isnanf(fn)) return x*fn; + if (!finitef(fn)) { + if(fn>(float)0.0) return x*fn; + else return x/(-fn); + } + if (rintf(fn)!=fn) return (fn-fn)/(fn-fn); +#if INT_MAX > 65000 + if ( fn > (float)65000.0) return scalbnf(x, 65000); + if (-fn > (float)65000.0) return scalbnf(x,-65000); +#else + if ( fn > (float)32000.0) return scalbnf(x, 32000); + if (-fn > (float)32000.0) return scalbnf(x,-32000); +#endif + return scalbnf(x,(int)fn); +#endif +} Index: s_ceil.c =================================================================== --- s_ceil.c (nonexistent) +++ s_ceil.c (revision 1765) @@ -0,0 +1,38 @@ + +/* @(#)z_ceil.c 1.0 98/08/13 */ +/***************************************************************** + * ceil + * + * Input: + * x - floating point value + * + * Output: + * Smallest integer greater than x. + * + * Description: + * This routine returns the smallest integer greater than x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +double +_DEFUN (ceil, (double), + double x) +{ + double f, y; + + y = modf (x, &f); + + if (y == 0.0) + return (x); + else if (x > -1.0 && x < 1.0) + return (x > 0 ? 1.0 : 0.0); + else + return (x > 0 ? f + 1.0 : f); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: sf_sineh.c =================================================================== --- sf_sineh.c (nonexistent) +++ sf_sineh.c (revision 1765) @@ -0,0 +1,110 @@ + +/* @(#)z_sinehf.c 1.0 98/08/13 */ +/****************************************************************** + * The following routines are coded directly from the algorithms + * and coefficients given in "Software Manual for the Elementary + * Functions" by William J. Cody, Jr. and William Waite, Prentice + * Hall, 1980. + ******************************************************************/ +/****************************************************************** + * Hyperbolic Sine + * + * Input: + * x - floating point value + * + * Output: + * hyperbolic sine of x + * + * Description: + * This routine calculates hyperbolic sines. + * + *****************************************************************/ + +#include +#include "fdlibm.h" +#include "zmath.h" + +static const float q[] = { -0.428277109e+2 }; +static const float p[] = { -0.713793159e+1, + -0.190333399 }; +static const float LNV = 0.6931610107; +static const float INV_V2 = 0.2499930850; +static const float V_OVER2_MINUS1 = 0.1383027787e-4; + +float +_DEFUN (sinehf, (float, int), + float x _AND + int cosineh) +{ + float y, f, P, Q, R, res, z, w; + int sgn = 1; + float WBAR = 18.55; + + /* Check for special values. */ + switch (numtestf (x)) + { + case NAN: + errno = EDOM; + return (x); + case INF: + errno = ERANGE; + return (ispos (x) ? z_infinity_f.f : -z_infinity_f.f); + } + + y = fabs (x); + + if (!cosineh && x < 0.0) + sgn = -1; + + if ((y > 1.0 && !cosineh) || cosineh) + { + if (y > BIGX) + { + w = y - LNV; + + /* Check for w > maximum here. */ + if (w > BIGX) + { + errno = ERANGE; + return (x); + } + + z = exp (w); + + if (w > WBAR) + res = z * (V_OVER2_MINUS1 + 1.0); + } + + else + { + z = exp (y); + if (cosineh) + res = (z + 1 / z) / 2.0; + else + res = (z - 1 / z) / 2.0; + } + + if (sgn < 0) + res = -res; + } + else + { + /* Check for y being too small. */ + if (y < z_rooteps_f) + { + res = x; + } + /* Calculate the Taylor series. */ + else + { + f = x * x; + Q = f + q[0]; + P = p[1] * f + p[0]; + R = f * (P / Q); + + res = x + x * R; + } + } + + return (res); +}

sf_sineh.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: s_log.c =================================================================== --- s_log.c (nonexistent) +++ s_log.c (revision 1765) @@ -0,0 +1,29 @@ + +/* @(#)z_log.c 1.0 98/08/13 */ +/****************************************************************** + * Logarithm + * + * Input: + * x - floating point value + * + * Output: + * natural logarithm of x + * + * Description: + * This routine returns the natural logarithm of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +double +_DEFUN (log, (double), + double x) +{ + return (logarithm (x, 0)); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: sf_fabs.c =================================================================== --- sf_fabs.c (nonexistent) +++ sf_fabs.c (revision 1765) @@ -0,0 +1,45 @@ + +/* @(#)z_fabsf.c 1.0 98/08/13 */ +/****************************************************************** + * Floating-Point Absolute Value + * + * Input: + * x - floating-point number + * + * Output: + * absolute value of x + * + * Description: + * fabs computes the absolute value of a floating point number. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +float +_DEFUN (fabsf, (float), + float x) +{ + switch (numtestf (x)) + { + case NAN: + errno = EDOM; + return (x); + case INF: + errno = ERANGE; + return (x); + case 0: + return (0.0); + default: + return (x < 0.0 ? -x : x); + } +} + +#ifdef _DOUBLE_IS_32BITS +double fabs (double x) +{ + return (double) fabsf ((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */

sf_fabs.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: s_tan.c =================================================================== --- s_tan.c (nonexistent) +++ s_tan.c (revision 1765) @@ -0,0 +1,139 @@ + +/* @(#)z_tan.c 1.0 98/08/13 */ +/****************************************************************** + * The following routines are coded directly from the algorithms + * and coefficients given in "Software Manual for the Elementary + * Functions" by William J. Cody, Jr. and William Waite, Prentice + * Hall, 1980. + ******************************************************************/ + +/* +FUNCTION + <>, <>---tangent + +INDEX +tan +INDEX +tanf + +ANSI_SYNOPSIS + #include + double tan(double <[x]>); + float tanf(float <[x]>); + +TRAD_SYNOPSIS + #include + double tan(<[x]>) + double <[x]>; + + float tanf(<[x]>) + float <[x]>; + + +DESCRIPTION +<> computes the tangent of the argument <[x]>. +Angles are specified in radians. + +<> is identical, save that it takes and returns <> values. + +RETURNS +The tangent of <[x]> is returned. + +PORTABILITY +<> is ANSI. <> is an extension. +*/ + +/****************************************************************** + * Tangent + * + * Input: + * x - floating point value + * + * Output: + * tangent of x + * + * Description: + * This routine calculates the tangent of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +static const double TWO_OVER_PI = 0.63661977236758134308; +static const double p[] = { -0.13338350006421960681, + 0.34248878235890589960e-2, + -0.17861707342254426711e-4 }; +static const double q[] = { -0.46671683339755294240, + 0.25663832289440112864e-1, + -0.31181531907010027307e-3, + 0.49819433993786512270e-6 }; + +double +_DEFUN (tan, (double), + double x) +{ + double y, f, g, XN, xnum, xden, res; + int N; + + /* Check for special values. */ + switch (numtest (x)) + { + case NAN: + errno = EDOM; + return (x); + case INF: + errno = EDOM; + return (z_notanum.d); + } + + y = fabs (x); + + /* Check for values that are out of our range. */ + if (y > 105414357.0) + { + errno = ERANGE; + return (y); + } + + if (x < 0.0) + N = (int) (x * TWO_OVER_PI - 0.5); + else + N = (int) (x * TWO_OVER_PI + 0.5); + + XN = (double) N; + + f = x - N * __PI_OVER_TWO; + + /* Check for values that are too small. */ + if (-z_rooteps < f && f < z_rooteps) + { + xnum = f; + xden = 1.0; + } + + /* Calculate the polynomial. */ + else + { + g = f * f; + + xnum = f * ((p[2] * g + p[1]) * g + p[0]) * g + f; + xden = (((q[3] * g + q[2]) * g + q[1]) * g + q[0]) * g + 1.0; + } + + if (N & 1) + { + xnum = -xnum; + res = xden / xnum; + } + else + { + res = xnum / xden; + } + + return (res); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: sf_ceil.c =================================================================== --- sf_ceil.c (nonexistent) +++ sf_ceil.c (revision 1765) @@ -0,0 +1,42 @@ + +/* @(#)z_ceilf.c 1.0 98/08/13 */ +/***************************************************************** + * ceil + * + * Input: + * x - floating point value + * + * Output: + * Smallest integer greater than x. + * + * Description: + * This routine returns the smallest integer greater than x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +float +_DEFUN (ceilf, (float), + float x) +{ + float f, y; + + y = modff (x, &f); + + if (y == 0.0) + return (x); + else if (x > -1.0 && x < 1.0) + return (x > 0 ? 1.0 : 0.0); + else + return (x > 0 ? f + 1.0 : f); +} + +#ifdef _DOUBLE_IS_32BITS +double ceil (double x) +{ + return (double) ceilf ((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */

sf_ceil.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: s_cos.c =================================================================== --- s_cos.c (nonexistent) +++ s_cos.c (revision 1765) @@ -0,0 +1,29 @@ + +/* @(#)z_cos.c 1.0 98/08/13 */ +/****************************************************************** + * Cosine + * + * Input: + * x - floating point value + * + * Output: + * cosine of x + * + * Description: + * This routine returns the cosine of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +double +_DEFUN (cos, (double), + double x) +{ + return (sine (x, 1)); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: sf_log10.c =================================================================== --- sf_log10.c (nonexistent) +++ sf_log10.c (revision 1765) @@ -0,0 +1,34 @@ + +/* @(#)z_log10f.c 1.0 98/08/13 */ +/****************************************************************** + * Logarithm + * + * Input: + * x - floating point value + * + * Output: + * logarithm of x + * + * Description: + * This routine returns the logarithm of x (base 10). + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +float +_DEFUN (log10f, (float), + float x) +{ + return (logarithmf (x, 1)); +} + +#ifdef _DOUBLE_IS_32BITS + +double log10 (double x) +{ + return (double) log10f ((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */

sf_log10.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: s_sqrt.c =================================================================== --- s_sqrt.c (nonexistent) +++ s_sqrt.c (revision 1765) @@ -0,0 +1,129 @@ + +/* @(#)z_sqrt.c 1.0 98/08/13 */ +/***************************************************************** + * The following routines are coded directly from the algorithms + * and coefficients given in "Software Manual for the Elementary + * Functions" by William J. Cody, Jr. and William Waite, Prentice + * Hall, 1980. + *****************************************************************/ + +/* +FUNCTION + <>, <>---positive square root + +INDEX + sqrt +INDEX + sqrtf + +ANSI_SYNOPSIS + #include + double sqrt(double <[x]>); + float sqrtf(float <[x]>); + +TRAD_SYNOPSIS + #include + double sqrt(<[x]>); + float sqrtf(<[x]>); + +DESCRIPTION + <> computes the positive square root of the argument. + +RETURNS + On success, the square root is returned. If <[x]> is real and + positive, then the result is positive. If <[x]> is real and + negative, the global value <> is set to <> (domain error). + + +PORTABILITY + <> is ANSI C. <> is an extension. +*/ + +/****************************************************************** + * Square Root + * + * Input: + * x - floating point value + * + * Output: + * square-root of x + * + * Description: + * This routine performs floating point square root. + * + * The initial approximation is computed as + * y0 = 0.41731 + 0.59016 * f + * where f is a fraction such that x = f * 2^exp. + * + * Three Newton iterations in the form of Heron's formula + * are then performed to obtain the final value: + * y[i] = (y[i-1] + f / y[i-1]) / 2, i = 1, 2, 3. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +double +_DEFUN (sqrt, (double), + double x) +{ + double f, y; + int exp, i, odd; + + /* Check for special values. */ + switch (numtest (x)) + { + case NAN: + errno = EDOM; + return (x); + case INF: + if (ispos (x)) + { + errno = EDOM; + return (z_notanum.d); + } + else + { + errno = ERANGE; + return (z_infinity.d); + } + } + + /* Initial checks are performed here. */ + if (x == 0.0) + return (0.0); + if (x < 0) + { + errno = EDOM; + return (z_notanum.d); + } + + /* Find the exponent and mantissa for the form x = f * 2^exp. */ + f = frexp (x, &exp); + + odd = exp & 1; + + /* Get the initial approximation. */ + y = 0.41731 + 0.59016 * f; + + f /= 2.0; + /* Calculate the remaining iterations. */ + for (i = 0; i < 3; ++i) + y = y / 2.0 + f / y; + + /* Calculate the final value. */ + if (odd) + { + y *= __SQRT_HALF; + exp++; + } + exp >>= 1; + y = ldexp (y, exp); + + return (y); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: s_sin.c =================================================================== --- s_sin.c (nonexistent) +++ s_sin.c (revision 1765) @@ -0,0 +1,29 @@ + +/* @(#)z_sin.c 1.0 98/08/13 */ +/****************************************************************** + * Sine + * + * Input: + * x - floating point value + * + * Output: + * sine of x + * + * Description: + * This routine returns the sine of x. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +double +_DEFUN (sin, (double), + double x) +{ + return (sine (x, 0)); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: s_tanh.c =================================================================== --- s_tanh.c (nonexistent) +++ s_tanh.c (revision 1765) @@ -0,0 +1,117 @@ + +/* @(#)z_tanh.c 1.0 98/08/13 */ +/***************************************************************** + * The following routines are coded directly from the algorithms + * and coefficients given in "Software Manual for the Elementary + * Functions" by William J. Cody, Jr. and William Waite, Prentice + * Hall, 1980. + *****************************************************************/ + +/* + +FUNCTION + <>, <>---hyperbolic tangent + +INDEX +tanh +INDEX +tanhf + +ANSI_SYNOPSIS + #include + double tanh(double <[x]>); + float tanhf(float <[x]>); + +TRAD_SYNOPSIS + #include + double tanh(<[x]>) + double <[x]>; + + float tanhf(<[x]>) + float <[x]>; + + +DESCRIPTION + +<> computes the hyperbolic tangent of +the argument <[x]>. Angles are specified in radians. + +<)>> is defined as +. sinh(<[x]>)/cosh(<[x]>) + +<> is identical, save that it takes and returns <> values. + +RETURNS +The hyperbolic tangent of <[x]> is returned. + +PORTABILITY +<> is ANSI C. <> is an extension. + +*/ + +/****************************************************************** + * Hyperbolic Tangent + * + * Input: + * x - floating point value + * + * Output: + * hyperbolic tangent of x + * + * Description: + * This routine calculates hyperbolic tangent. + * + *****************************************************************/ + +#include +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +static const double LN3_OVER2 = 0.54930614433405484570; +static const double p[] = { -0.16134119023996228053e+4, + -0.99225929672236083313e+2, + -0.96437492777225469787 }; +static const double q[] = { 0.48402357071988688686e+4, + 0.22337720718962312926e+4, + 0.11274474380534949335e+3 }; + +double +_DEFUN (tanh, (double), + double x) +{ + double f, res, g, P, Q, R; + + f = fabs (x); + + /* Check if the input is too big. */ + if (f > BIGX) + res = 1.0; + + else if (f > LN3_OVER2) + res = 1.0 - 2.0 / (exp (2 * f) + 1.0); + + /* Check if the input is too small. */ + else if (f < z_rooteps) + res = f; + + /* Calculate the Taylor series. */ + else + { + g = f * f; + + P = (p[2] * g + p[1]) * g + p[0]; + Q = ((g + q[2]) * g + q[1]) * g + q[0]; + R = g * (P / Q); + + res = f + f * R; + } + + if (x < 0.0) + res = -res; + + return (res); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: sf_frexp.c =================================================================== --- sf_frexp.c (nonexistent) +++ sf_frexp.c (revision 1765) @@ -0,0 +1,58 @@ + +/* @(#)z_frexpf.c 1.0 98/08/13 */ +/****************************************************************** + * frexp + * + * Input: + * d - floating point value + * exp - exponent value + * + * Output: + * A floating point value in the range [0.5, 1). + * + * Description: + * This routine breaks a floating point value into a number f and + * an exponent exp such that d = f * 2 ^ exp. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +float frexpf (float d, int *exp) +{ + float f; + __int32_t wf, wd; + + GET_FLOAT_WORD (wd, d); + + /* Get the exponent. */ + *exp = ((wd & 0x7f800000) >> 23) - 126; + + /* Get the mantissa. */ + wf = wd & 0x7fffff; + wf |= 0x3f000000; + + SET_FLOAT_WORD (f, wf); + + /* Check for special values. */ + switch (numtestf (f)) + { + case NAN: + case INF: + errno = EDOM; + *exp = 0; + return (f); + } + + return (f); +} + +#ifdef _DOUBLE_IS_32BITS + +double frexp (double x, int *exp) +{ + return (double) frexpf ((float) x, exp); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */

sf_frexp.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: sf_logarithm.c =================================================================== --- sf_logarithm.c (nonexistent) +++ sf_logarithm.c (revision 1765) @@ -0,0 +1,72 @@ + +/* @(#)z_logarithmf.c 1.0 98/08/13 */ +/****************************************************************** + * The following routines are coded directly from the algorithms + * and coefficients given in "Software Manual for the Elementary + * Functions" by William J. Cody, Jr. and William Waite, Prentice + * Hall, 1980. + ******************************************************************/ +/****************************************************************** + * Logarithm + * + * Input: + * x - floating point value + * ten - indicates base ten numbers + * + * Output: + * logarithm of x + * + * Description: + * This routine calculates logarithms. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +static const float a[] = { -0.5527074855 }; +static const float b[] = { -0.6632718214e+1 }; +static const float C1 = 0.693145752; +static const float C2 = 1.428606820e-06; +static const float C3 = 0.4342944819; + +float +_DEFUN (logarithmf, (float, int), + float x _AND + int ten) +{ + int N; + float f, w, z; + + /* Check for domain error here. */ + if (x <= 0.0) + { + errno = ERANGE; + return (z_notanum_f.f); + } + + /* Get the exponent and mantissa where x = f * 2^N. */ + f = frexpf (x, &N); + + z = f - 0.5; + + if (f > __SQRT_HALF) + z = (z - 0.5) / (f * 0.5 + 0.5); + else + { + N--; + z /= (z * 0.5 + 0.5); + } + w = z * z; + + /* Use Newton's method with 4 terms. */ + z += z * w * (a[0]) / ((w + 1.0) * w + b[0]); + + if (N != 0) + z = (N * C2 + z) + N * C1; + + if (ten) + z *= C3; + + return (z); +}

sf_logarithm.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: s_ispos.c =================================================================== --- s_ispos.c (nonexistent) +++ s_ispos.c (revision 1765) @@ -0,0 +1,35 @@ + +/* @(#)z_ispos.c 1.0 98/08/13 */ +/****************************************************************** + * Numtest + * + * Input: + * x - pointer to a floating point value + * + * Output: + * An integer that indicates if the number is positive. + * + * Description: + * This routine returns an integer that indicates if the number + * passed in is positive (1) or negative (0). + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +int ispos (double x) +{ + __uint32_t hx; + + GET_HIGH_WORD (hx, x); + + if (hx & 0x80000000) + return (0); + else + return (1); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: s_asine.c =================================================================== --- s_asine.c (nonexistent) +++ s_asine.c (revision 1765) @@ -0,0 +1,186 @@ + +/* @(#)z_asine.c 1.0 98/08/13 */ +/****************************************************************** + * The following routines are coded directly from the algorithms + * and coefficients given in "Software Manual for the Elementary + * Functions" by William J. Cody, Jr. and William Waite, Prentice + * Hall, 1980. + ******************************************************************/ + +/* +FUNCTION + <>, <>, <>, <>, <>, <>---arc sine or cosine + +INDEX + asin +INDEX + asinf +INDEX + acos +INDEX + acosf +INDEX + asine +INDEX + asinef + +ANSI_SYNOPSIS + #include + double asine(double <[x]>); + float asinef(float <[x]>); + double asin(double <[x]>); + float asinf(float <[x]>); + double acos(double <[x]>); + float acosf(float <[x]>); + +TRAD_SYNOPSIS + #include + double asine(<[x]>); + double <[x]>; + + float asinef(<[x]>); + float <[x]>; + + double asin(<[x]>) + double <[x]>; + + float asinf(<[x]>) + float <[x]>; + + double acos(<[x]>) + double <[x]>; + + float acosf(<[x]>) + float <[x]>; + +DESCRIPTION + +<> computes the inverse sine or cosine of the argument <[x]>. +Arguments to <> and <> must be in the range @minus{}1 to 1. + +<> and <> are identical to <> and <>, other +than taking and returning floats. + +RETURNS +@ifinfo +<> and <> return values in radians, in the range of -pi/2 to pi/2. +@end ifinfo +@tex +<> and <> return values in radians, in the range of $-\pi/2$ to $\pi/2$. +@end tex + +If <[x]> is not in the range @minus{}1 to 1, <> and <> +return NaN (not a number), set the global variable <> to +<>, and issue a <> message. + +*/ + +/****************************************************************** + * Arcsine + * + * Input: + * x - floating point value + * acosine - indicates acos calculation + * + * Output: + * Arcsine of x. + * + * Description: + * This routine calculates arcsine / arccosine. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +static const double p[] = { -0.27368494524164255994e+2, + 0.57208227877891731407e+2, + -0.39688862997404877339e+2, + 0.10152522233806463645e+2, + -0.69674573447350646411 }; +static const double q[] = { -0.16421096714498560795e+3, + 0.41714430248260412556e+3, + -0.38186303361750149284e+3, + 0.15095270841030604719e+3, + -0.23823859153670238830e+2 }; +static const double a[] = { 0.0, 0.78539816339744830962 }; +static const double b[] = { 1.57079632679489661923, 0.78539816339744830962 }; + +double +_DEFUN (asine, (double, int), + double x _AND + int acosine) +{ + int flag, i; + int branch = 0; + double g, res, R, P, Q, y; + + /* Check for special values. */ + i = numtest (x); + if (i == NAN || i == INF) + { + errno = EDOM; + if (i == NAN) + return (x); + else + return (z_infinity.d); + } + + y = fabs (x); + flag = acosine; + + if (y > 0.5) + { + i = 1 - flag; + + /* Check for range error. */ + if (y > 1.0) + { + errno = ERANGE; + return (z_notanum.d); + } + + g = (1 - y) / 2.0; + y = -2 * sqrt (g); + branch = 1; + } + else + { + i = flag; + if (y < z_rooteps) + res = y; + else + g = y * y; + } + + if (y >= z_rooteps || branch == 1) + { + /* Calculate the Taylor series. */ + P = ((((p[4] * g + p[3]) * g + p[2]) * g + p[1]) * g + p[0]) * g; + Q = ((((g + q[4]) * g + q[3]) * g + q[2]) * g + q[1]) * g + q[0]; + R = P / Q; + + res = y + y * R; + } + + /* Calculate asine or acose. */ + if (flag == 0) + { + res = (a[i] + res) + a[i]; + if (x < 0.0) + res = -res; + } + else + { + if (x < 0.0) + res = (b[i] + res) + b[i]; + else + res = (a[i] - res) + a[i]; + } + + return (res); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: sf_sqrt.c =================================================================== --- sf_sqrt.c (nonexistent) +++ sf_sqrt.c (revision 1765) @@ -0,0 +1,100 @@ + +/* @(#)z_sqrtf.c 1.0 98/08/13 */ +/***************************************************************** + * The following routines are coded directly from the algorithms + * and coefficients given in "Software Manual for the Elementary + * Functions" by William J. Cody, Jr. and William Waite, Prentice + * Hall, 1980. + *****************************************************************/ +/****************************************************************** + * Square Root + * + * Input: + * x - floating point value + * + * Output: + * square-root of x + * + * Description: + * This routine performs floating point square root. + * + * The initial approximation is computed as + * y0 = 0.41731 + 0.59016 * f + * where f is a fraction such that x = f * 2^exp. + * + * Three Newton iterations in the form of Heron's formula + * are then performed to obtain the final value: + * y[i] = (y[i-1] + f / y[i-1]) / 2, i = 1, 2, 3. + * + *****************************************************************/ + +#include "fdlibm.h" +#include "zmath.h" + +float +_DEFUN (sqrtf, (float), + float x) +{ + float f, y; + int exp, i, odd; + + /* Check for special values. */ + switch (numtestf (x)) + { + case NAN: + errno = EDOM; + return (x); + case INF: + if (isposf (x)) + { + errno = EDOM; + return (z_notanum_f.f); + } + else + { + errno = ERANGE; + return (z_infinity_f.f); + } + } + + /* Initial checks are performed here. */ + if (x == 0.0) + return (0.0); + if (x < 0) + { + errno = EDOM; + return (z_notanum_f.f); + } + + /* Find the exponent and mantissa for the form x = f * 2^exp. */ + f = frexpf (x, &exp); + odd = exp & 1; + + /* Get the initial approximation. */ + y = 0.41731 + 0.59016 * f; + + f *= 0.5; + /* Calculate the remaining iterations. */ + for (i = 0; i < 2; ++i) + y = y * 0.5 + f / y; + + /* Calculate the final value. */ + if (odd) + { + y *= __SQRT_HALF; + exp++; + } + exp >>= 1; + y = ldexpf (y, exp); + + return (y); +} + +#ifdef _DOUBLE_IS_32BITS + +double sqrt (double x) +{ + return (double) sqrtf ((float) x); +} + +#endif /* _DOUBLE_IS_32BITS */

sf_sqrt.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: sf_tanh.c =================================================================== --- sf_tanh.c (nonexistent) +++ sf_tanh.c (revision 1765) @@ -0,0 +1,77 @@ + +/* @(#)z_tanhf.c 1.0 98/08/13 */ +/***************************************************************** + * The following routines are coded directly from the algorithms + * and coefficients given in "Software Manual for the Elementary + * Functions" by William J. Cody, Jr. and William Waite, Prentice + * Hall, 1980. + *****************************************************************/ +/****************************************************************** + * Hyperbolic Tangent + * + * Input: + * x - floating point value + * + * Output: + * hyperbolic tangent of x + * + * Description: + * This routine calculates hyperbolic tangent. + * + *****************************************************************/ + +#include +#include "fdlibm.h" +#include "zmath.h" + +static const float LN3_OVER2 = 0.5493061443; +static const float p[] = { -0.2059432032, + -0.0009577527 }; +static const float q[] = { 0.6178299136, + 0.25 }; + +float +_DEFUN (tanhf, (float), + float x) +{ + float f, res, g, P, Q, R; + + f = fabsf (x); + + /* Check if the input is too big. */ + if (f > BIGX) + res = 1.0; + + else if (f > LN3_OVER2) + res = 1.0 - 2.0 / (exp (2 * f) + 1.0); + + /* Check if the input is too small. */ + else if (f < z_rooteps_f) + res = f; + + /* Calculate the Taylor series. */ + else + { + g = f * f; + + P = p[1] * g + p[0]; + Q = (g + q[1]) * g + q[0]; + R = g * (P / Q); + + res = f + f * R; + } + + if (x < 0.0) + res = -res; + + return (res); +} + +#ifdef _DOUBLE_IS_32BITS + +double tanh (double x) +{ + return (double) tanhf ((float) x); +} + +#endif _DOUBLE_IS_32BITS

sf_tanh.c Property changes : Added: svn:executable ## -0,0 +1 ## +* \ No newline at end of property Index: s_atangent.c =================================================================== --- s_atangent.c (nonexistent) +++ s_atangent.c (revision 1765) @@ -0,0 +1,213 @@ + +/* @(#)z_atangent.c 1.0 98/08/13 */ +/****************************************************************** + * The following routines are coded directly from the algorithms + * and coefficients given in "Software Manual for the Elementary + * Functions" by William J. Cody, Jr. and William Waite, Prentice + * Hall, 1980. + ******************************************************************/ + +/* +FUNCTION + <>, <>, <>, <>, <>, <>---arc tangent + +INDEX + atan2 +INDEX + atan2f +INDEX + atan +INDEX + atanf + +ANSI_SYNOPSIS + #include + double atan(double <[x]>); + float atan(float <[x]>); + double atan2(double <[y]>,double <[x]>); + float atan2f(float <[y]>,float <[x]>); + +TRAD_SYNOPSIS + #include + double atan2(<[y]>,<[x]>); + double <[y]>; + double <[x]>; + + float atan2f(<[y]>,<[x]>); + float <[y]>; + float <[x]>; + + #include + double atan(<[x]>); + double <[x]>; + + float atanf(<[x]>); + float <[x]>; + +DESCRIPTION + +<> computes the inverse tangent (arc tangent) of y / x. + +<> is identical to <>, save that it operates on <>. + +<> computes the inverse tangent (arc tangent) of the input value. + +<> is identical to <>, save that it operates on <>. + +RETURNS +@ifinfo +<> returns a value in radians, in the range of -pi/2 to pi/2. +<> returns a value in radians, in the range of -pi/2 to pi/2. +@end ifinfo +@tex +<> returns a value in radians, in the range of $-\pi/2$ to $\pi/2$. +<> returns a value in radians, in the range of $-\pi/2$ to $\pi/2$. +@end tex + +PORTABILITY +<> is ANSI C. <> is an extension. +<> is ANSI C. <> is an extension. + +*/ + +/****************************************************************** + * Arctangent + * + * Input: + * x - floating point value + * + * Output: + * arctangent of x + * + * Description: + * This routine calculates arctangents. + * + *****************************************************************/ +#include +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +static const double ROOT3 = 1.73205080756887729353; +static const double a[] = { 0.0, 0.52359877559829887308, 1.57079632679489661923, + 1.04719755119659774615 }; +static const double q[] = { 0.41066306682575781263e+2, + 0.86157349597130242515e+2, + 0.59578436142597344465e+2, + 0.15024001160028576121e+2 }; +static const double p[] = { -0.13688768894191926929e+2, + -0.20505855195861651981e+2, + -0.84946240351320683534e+1, + -0.83758299368150059274 }; + +double +_DEFUN (atangent, (double, double, double, int), + double x _AND + double v _AND + double u _AND + int arctan2) +{ + double f, g, R, P, Q, A, res; + int N; + int branch = 0; + int expv, expu; + + /* Preparation for calculating arctan2. */ + if (arctan2) + { + if (u == 0.0) + if (v == 0.0) + { + errno = ERANGE; + return (z_notanum.d); + } + else + { + branch = 1; + res = __PI_OVER_TWO; + } + + if (!branch) + { + int e; + /* Get the exponent values of the inputs. */ + g = frexp (v, &expv); + g = frexp (u, &expu); + + /* See if a divide will overflow. */ + e = expv - expu; + if (e > DBL_MAX_EXP) + { + branch = 1; + res = __PI_OVER_TWO; + } + + /* Also check for underflow. */ + else if (e < DBL_MIN_EXP) + { + branch = 2; + res = 0.0; + } + } + } + + if (!branch) + { + if (arctan2) + f = fabs (v / u); + else + f = fabs (x); + + if (f > 1.0) + { + f = 1.0 / f; + N = 2; + } + else + N = 0; + + if (f > (2.0 - ROOT3)) + { + A = ROOT3 - 1.0; + f = (((A * f - 0.5) - 0.5) + f) / (ROOT3 + f); + N++; + } + + /* Check for values that are too small. */ + if (-z_rooteps < f && f < z_rooteps) + res = f; + + /* Calculate the Taylor series. */ + else + { + g = f * f; + P = (((p[3] * g + p[2]) * g + p[1]) * g + p[0]) * g; + Q = (((g + q[3]) * g + q[2]) * g + q[1]) * g + q[0]; + R = P / Q; + + res = f + f * R; + } + + if (N > 1) + res = -res; + + res += a[N]; + } + + if (arctan2) + { + if (u < 0.0 || branch == 2) + res = __PI - res; + if (v < 0.0 || branch == 1) + res = -res; + } + else if (x < 0.0) + { + res = -res; + } + + return (res); +} + +#endif /* _DOUBLE_IS_32BITS */ Index: er_lgamma.c =================================================================== --- er_lgamma.c (nonexistent) +++ er_lgamma.c (revision 1765) @@ -0,0 +1,422 @@ + +/* @(#)er_lgamma.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + */ + +/* +FUNCTION + <>, <>, <>, <>, <>, + <>, <>, <>---logarithmic gamma + function +INDEX +gamma +INDEX +gammaf +INDEX +lgamma +INDEX +lgammaf +INDEX +gamma_r +INDEX +gammaf_r +INDEX +lgamma_r +INDEX +lgammaf_r + +ANSI_SYNOPSIS +#include +double gamma(double <[x]>); +float gammaf(float <[x]>); +double lgamma(double <[x]>); +float lgammaf(float <[x]>); +double gamma_r(double <[x]>, int *<[signgamp]>); +float gammaf_r(float <[x]>, int *<[signgamp]>); +double lgamma_r(double <[x]>, int *<[signgamp]>); +float lgammaf_r(float <[x]>, int *<[signgamp]>); + +TRAD_SYNOPSIS +#include +double gamma(<[x]>) +double <[x]>; +float gammaf(<[x]>) +float <[x]>; +double lgamma(<[x]>) +double <[x]>; +float lgammaf(<[x]>) +float <[x]>; +double gamma_r(<[x]>, <[signgamp]>) +double <[x]>; +int <[signgamp]>; +float gammaf_r(<[x]>, <[signgamp]>) +float <[x]>; +int <[signgamp]>; +double lgamma_r(<[x]>, <[signgamp]>) +double <[x]>; +int <[signgamp]>; +float lgammaf_r(<[x]>, <[signgamp]>) +float <[x]>; +int <[signgamp]>; + +DESCRIPTION +<> calculates +@tex +$\mit ln\bigl(\Gamma(x)\bigr)$, +@end tex +the natural logarithm of the gamma function of <[x]>. The gamma function +(<))>>) is a generalization of factorial, and retains +the property that +@ifinfo +<> is equivalent to <>. +@end ifinfo +@tex +$\mit \Gamma(N)\equiv N\times\Gamma(N-1)$. +@end tex +Accordingly, the results of the gamma function itself grow very +quickly. <> is defined as +@tex +$\mit ln\bigl(\Gamma(x)\bigr)$ rather than simply $\mit \Gamma(x)$ +@end tex +@ifinfo +the natural log of the gamma function, rather than the gamma function +itself, +@end ifinfo +to extend the useful range of results representable. + +The sign of the result is returned in the global variable <>, +which is declared in math.h. + +<> performs the same calculation as <>, but uses and +returns <> values. + +<> and <> are alternate names for <> and +<>. The use of <> instead of <> is a reminder +that these functions compute the log of the gamma function, rather +than the gamma function itself. + +The functions <>, <>, <>, and +<> are just like <>, <>, <>, and +<>, respectively, but take an additional argument. This +additional argument is a pointer to an integer. This additional +argument is used to return the sign of the result, and the global +variable <> is not used. These functions may be used for +reentrant calls (but they will still set the global variable <> +if an error occurs). + +RETURNS +Normally, the computed result is returned. + +When <[x]> is a nonpositive integer, <> returns <> +and <> is set to <>. If the result overflows, <> +returns <> and <> is set to <>. + +You can modify this error treatment using <>. + +PORTABILITY +Neither <> nor <> is ANSI C. */ + +/* lgamma_r(x, signgamp) + * Reentrant version of the logarithm of the Gamma function + * with user provide pointer for the sign of Gamma(x). + * + * Method: + * 1. Argument Reduction for 0 < x <= 8 + * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may + * reduce x to a number in [1.5,2.5] by + * lgamma(1+s) = log(s) + lgamma(s) + * for example, + * lgamma(7.3) = log(6.3) + lgamma(6.3) + * = log(6.3*5.3) + lgamma(5.3) + * = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3) + * 2. Polynomial approximation of lgamma around its + * minimun ymin=1.461632144968362245 to maintain monotonicity. + * On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use + * Let z = x-ymin; + * lgamma(x) = -1.214862905358496078218 + z^2*poly(z) + * where + * poly(z) is a 14 degree polynomial. + * 2. Rational approximation in the primary interval [2,3] + * We use the following approximation: + * s = x-2.0; + * lgamma(x) = 0.5*s + s*P(s)/Q(s) + * with accuracy + * |P/Q - (lgamma(x)-0.5s)| < 2**-61.71 + * Our algorithms are based on the following observation + * + * zeta(2)-1 2 zeta(3)-1 3 + * lgamma(2+s) = s*(1-Euler) + --------- * s - --------- * s + ... + * 2 3 + * + * where Euler = 0.5771... is the Euler constant, which is very + * close to 0.5. + * + * 3. For x>=8, we have + * lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+.... + * (better formula: + * lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...) + * Let z = 1/x, then we approximation + * f(z) = lgamma(x) - (x-0.5)(log(x)-1) + * by + * 3 5 11 + * w = w0 + w1*z + w2*z + w3*z + ... + w6*z + * where + * |w - f(z)| < 2**-58.74 + * + * 4. For negative x, since (G is gamma function) + * -x*G(-x)*G(x) = pi/sin(pi*x), + * we have + * G(x) = pi/(sin(pi*x)*(-x)*G(-x)) + * since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0 + * Hence, for x<0, signgam = sign(sin(pi*x)) and + * lgamma(x) = log(|Gamma(x)|) + * = log(pi/(|x*sin(pi*x)|)) - lgamma(-x); + * Note: one should avoid compute pi*(-x) directly in the + * computation of sin(pi*(-x)). + * + * 5. Special Cases + * lgamma(2+s) ~ s*(1-Euler) for tiny s + * lgamma(1)=lgamma(2)=0 + * lgamma(x) ~ -log(x) for tiny x + * lgamma(0) = lgamma(inf) = inf + * lgamma(-integer) = +-inf + * + */ + +#include "fdlibm.h" + +#ifdef __STDC__ +static const double +#else +static double +#endif +two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ +half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ +one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ +pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ +a0 = 7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */ +a1 = 3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */ +a2 = 6.73523010531292681824e-02, /* 0x3FB13E00, 0x1A5562A7 */ +a3 = 2.05808084325167332806e-02, /* 0x3F951322, 0xAC92547B */ +a4 = 7.38555086081402883957e-03, /* 0x3F7E404F, 0xB68FEFE8 */ +a5 = 2.89051383673415629091e-03, /* 0x3F67ADD8, 0xCCB7926B */ +a6 = 1.19270763183362067845e-03, /* 0x3F538A94, 0x116F3F5D */ +a7 = 5.10069792153511336608e-04, /* 0x3F40B6C6, 0x89B99C00 */ +a8 = 2.20862790713908385557e-04, /* 0x3F2CF2EC, 0xED10E54D */ +a9 = 1.08011567247583939954e-04, /* 0x3F1C5088, 0x987DFB07 */ +a10 = 2.52144565451257326939e-05, /* 0x3EFA7074, 0x428CFA52 */ +a11 = 4.48640949618915160150e-05, /* 0x3F07858E, 0x90A45837 */ +tc = 1.46163214496836224576e+00, /* 0x3FF762D8, 0x6356BE3F */ +tf = -1.21486290535849611461e-01, /* 0xBFBF19B9, 0xBCC38A42 */ +/* tt = -(tail of tf) */ +tt = -3.63867699703950536541e-18, /* 0xBC50C7CA, 0xA48A971F */ +t0 = 4.83836122723810047042e-01, /* 0x3FDEF72B, 0xC8EE38A2 */ +t1 = -1.47587722994593911752e-01, /* 0xBFC2E427, 0x8DC6C509 */ +t2 = 6.46249402391333854778e-02, /* 0x3FB08B42, 0x94D5419B */ +t3 = -3.27885410759859649565e-02, /* 0xBFA0C9A8, 0xDF35B713 */ +t4 = 1.79706750811820387126e-02, /* 0x3F9266E7, 0x970AF9EC */ +t5 = -1.03142241298341437450e-02, /* 0xBF851F9F, 0xBA91EC6A */ +t6 = 6.10053870246291332635e-03, /* 0x3F78FCE0, 0xE370E344 */ +t7 = -3.68452016781138256760e-03, /* 0xBF6E2EFF, 0xB3E914D7 */ +t8 = 2.25964780900612472250e-03, /* 0x3F6282D3, 0x2E15C915 */ +t9 = -1.40346469989232843813e-03, /* 0xBF56FE8E, 0xBF2D1AF1 */ +t10 = 8.81081882437654011382e-04, /* 0x3F4CDF0C, 0xEF61A8E9 */ +t11 = -5.38595305356740546715e-04, /* 0xBF41A610, 0x9C73E0EC */ +t12 = 3.15632070903625950361e-04, /* 0x3F34AF6D, 0x6C0EBBF7 */ +t13 = -3.12754168375120860518e-04, /* 0xBF347F24, 0xECC38C38 */ +t14 = 3.35529192635519073543e-04, /* 0x3F35FD3E, 0xE8C2D3F4 */ +u0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */ +u1 = 6.32827064025093366517e-01, /* 0x3FE4401E, 0x8B005DFF */ +u2 = 1.45492250137234768737e+00, /* 0x3FF7475C, 0xD119BD6F */ +u3 = 9.77717527963372745603e-01, /* 0x3FEF4976, 0x44EA8450 */ +u4 = 2.28963728064692451092e-01, /* 0x3FCD4EAE, 0xF6010924 */ +u5 = 1.33810918536787660377e-02, /* 0x3F8B678B, 0xBF2BAB09 */ +v1 = 2.45597793713041134822e+00, /* 0x4003A5D7, 0xC2BD619C */ +v2 = 2.12848976379893395361e+00, /* 0x40010725, 0xA42B18F5 */ +v3 = 7.69285150456672783825e-01, /* 0x3FE89DFB, 0xE45050AF */ +v4 = 1.04222645593369134254e-01, /* 0x3FBAAE55, 0xD6537C88 */ +v5 = 3.21709242282423911810e-03, /* 0x3F6A5ABB, 0x57D0CF61 */ +s0 = -7.72156649015328655494e-02, /* 0xBFB3C467, 0xE37DB0C8 */ +s1 = 2.14982415960608852501e-01, /* 0x3FCB848B, 0x36E20878 */ +s2 = 3.25778796408930981787e-01, /* 0x3FD4D98F, 0x4F139F59 */ +s3 = 1.46350472652464452805e-01, /* 0x3FC2BB9C, 0xBEE5F2F7 */ +s4 = 2.66422703033638609560e-02, /* 0x3F9B481C, 0x7E939961 */ +s5 = 1.84028451407337715652e-03, /* 0x3F5E26B6, 0x7368F239 */ +s6 = 3.19475326584100867617e-05, /* 0x3F00BFEC, 0xDD17E945 */ +r1 = 1.39200533467621045958e+00, /* 0x3FF645A7, 0x62C4AB74 */ +r2 = 7.21935547567138069525e-01, /* 0x3FE71A18, 0x93D3DCDC */ +r3 = 1.71933865632803078993e-01, /* 0x3FC601ED, 0xCCFBDF27 */ +r4 = 1.86459191715652901344e-02, /* 0x3F9317EA, 0x742ED475 */ +r5 = 7.77942496381893596434e-04, /* 0x3F497DDA, 0xCA41A95B */ +r6 = 7.32668430744625636189e-06, /* 0x3EDEBAF7, 0xA5B38140 */ +w0 = 4.18938533204672725052e-01, /* 0x3FDACFE3, 0x90C97D69 */ +w1 = 8.33333333333329678849e-02, /* 0x3FB55555, 0x5555553B */ +w2 = -2.77777777728775536470e-03, /* 0xBF66C16C, 0x16B02E5C */ +w3 = 7.93650558643019558500e-04, /* 0x3F4A019F, 0x98CF38B6 */ +w4 = -5.95187557450339963135e-04, /* 0xBF4380CB, 0x8C0FE741 */ +w5 = 8.36339918996282139126e-04, /* 0x3F4B67BA, 0x4CDAD5D1 */ +w6 = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */ + +#ifdef __STDC__ +static const double zero= 0.00000000000000000000e+00; +#else +static double zero= 0.00000000000000000000e+00; +#endif + +#ifdef __STDC__ + static double sin_pi(double x) +#else + static double sin_pi(x) + double x; +#endif +{ + double y,z; + __int32_t n,ix; + + GET_HIGH_WORD(ix,x); + ix &= 0x7fffffff; + + if(ix<0x3fd00000) return __kernel_sin(pi*x,zero,0); + y = -x; /* x is assume negative */ + + /* + * argument reduction, make sure inexact flag not raised if input + * is an integer + */ + z = floor(y); + if(z!=y) { /* inexact anyway */ + y *= 0.5; + y = 2.0*(y - floor(y)); /* y = |x| mod 2.0 */ + n = (__int32_t) (y*4.0); + } else { + if(ix>=0x43400000) { + y = zero; n = 0; /* y must be even */ + } else { + if(ix<0x43300000) z = y+two52; /* exact */ + GET_LOW_WORD(n,z); + n &= 1; + y = n; + n<<= 2; + } + } + switch (n) { + case 0: y = __kernel_sin(pi*y,zero,0); break; + case 1: + case 2: y = __kernel_cos(pi*(0.5-y),zero); break; + case 3: + case 4: y = __kernel_sin(pi*(one-y),zero,0); break; + case 5: + case 6: y = -__kernel_cos(pi*(y-1.5),zero); break; + default: y = __kernel_sin(pi*(y-2.0),zero,0); break; + } + return -y; +} + + +#ifdef __STDC__ + double lgamma_r(double x, int *signgamp) +#else + double lgamma_r(x,signgamp) + double x; int *signgamp; +#endif +{ + double t,y,z,nadj,p,p1,p2,p3,q,r,w; + __int32_t i,hx,lx,ix; + + EXTRACT_WORDS(hx,lx,x); + + /* purge off +-inf, NaN, +-0, and negative arguments */ + *signgamp = 1; + ix = hx&0x7fffffff; + if(ix>=0x7ff00000) return x*x; + if((ix|lx)==0) return one/zero; + if(ix<0x3b900000) { /* |x|<2**-70, return -log(|x|) */ + if(hx<0) { + *signgamp = -1; + return -log(-x); + } else return -log(x); + } + if(hx<0) { + if(ix>=0x43300000) /* |x|>=2**52, must be -integer */ + return one/zero; + t = sin_pi(x); + if(t==zero) return one/zero; /* -integer */ + nadj = log(pi/fabs(t*x)); + if(t=0x3FE76944) {y = one-x; i= 0;} + else if(ix>=0x3FCDA661) {y= x-(tc-one); i=1;} + else {y = x; i=2;} + } else { + r = zero; + if(ix>=0x3FFBB4C3) {y=2.0-x;i=0;} /* [1.7316,2] */ + else if(ix>=0x3FF3B4C4) {y=x-tc;i=1;} /* [1.23,1.73] */ + else {y=x-one;i=2;} + } + switch(i) { + case 0: + z = y*y; + p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10)))); + p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11))))); + p = y*p1+p2; + r += (p-0.5*y); break; + case 1: + z = y*y; + w = z*y; + p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */ + p2 = t1+w*(t4+w*(t7+w*(t10+w*t13))); + p3 = t2+w*(t5+w*(t8+w*(t11+w*t14))); + p = z*p1-(tt-w*(p2+y*p3)); + r += (tf + p); break; + case 2: + p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5))))); + p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5)))); + r += (-0.5*y + p1/p2); + } + } + else if(ix<0x40200000) { /* x < 8.0 */ + i = (__int32_t)x; + t = zero; + y = x-(double)i; + p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); + q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); + r = half*y+p/q; + z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ + switch(i) { + case 7: z *= (y+6.0); /* FALLTHRU */ + case 6: z *= (y+5.0); /* FALLTHRU */ + case 5: z *= (y+4.0); /* FALLTHRU */ + case 4: z *= (y+3.0); /* FALLTHRU */ + case 3: z *= (y+2.0); /* FALLTHRU */ + r += log(z); break; + } + /* 8.0 <= x < 2**58 */ + } else if (ix < 0x43900000) { + t = log(x); + z = one/x; + y = z*z; + w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); + r = (x-half)*(t-one)+w; + } else + /* 2**58 <= x <= inf */ + r = x*(log(x)-one); + if(hx<0) r = nadj - r; + return r; +} Index: s_pow.c =================================================================== --- s_pow.c (nonexistent) +++ s_pow.c (revision 1765) @@ -0,0 +1,146 @@ + +/* @(#)z_pow.c 1.0 98/08/13 */ + +/* +FUNCTION + <>, <>---x to the power y +INDEX + pow +INDEX + powf + + +ANSI_SYNOPSIS + #include + double pow(double <[x]>, double <[y]>); + float pow(float <[x]>, float <[y]>); + +TRAD_SYNOPSIS + #include + double pow(<[x]>, <[y]>); + double <[x]>, <[y]>; + + float pow(<[x]>, <[y]>); + float <[x]>, <[y]>; + +DESCRIPTION + <> and <> calculate <[x]> raised to the exp1.0nt <[y]>. + @tex + (That is, $x^y$.) + @end tex + +RETURNS + On success, <> and <> return the value calculated. + + When the argument values would produce overflow, <> + returns <> and set <> to <>. If the + argument <[x]> passed to <> or <> is a negative + noninteger, and <[y]> is also not an integer, then <> + is set to <>. If <[x]> and <[y]> are both 0, then + <> and <> return <<1>>. + + You can modify error handling for these functions using <>. + +PORTABILITY + <> is ANSI C. <> is an extension. */ + +#include +#include "fdlibm.h" +#include "zmath.h" + +#ifndef _DOUBLE_IS_32BITS + +double pow (double x, double y) +{ + double d, t, r = 1.0; + int n, k, sign = 0; + __uint32_t px; + + GET_HIGH_WORD (px, x); + + k = modf (y, &d); + if (k == 0.0) + { + if (modf (ldexp (y, -1), &t)) + sign = 0; + else + sign = 1; + } + + if (x == 0.0 && y <= 0.0) + errno = EDOM; + + else if ((t = y * log (fabs (x))) >= BIGX) + { + errno = ERANGE; + if (px & 0x80000000) + { + if (!k) + { + errno = EDOM; + x = 0.0; + } + else if (sign) + x = -z_infinity.d; + else + x = z_infinity.d; + } + + else + x = z_infinity.d; + } + + else if (t < SMALLX) + { + errno = ERANGE; + x = 0.0; + } + + else + { + if ( k && fabs(d) <= 32767 ) + { + n = (int) d; + + if (sign = (n < 0)) + n = -n; + + while ( n > 0 ) + { + if ((unsigned int) n % 2) + r *= x; + x *= x; + n = (unsigned int) n / 2; + } + + if (sign) + r = 1.0 / r; + + return r; + } + + else + { + if ( px & 0x80000000 ) + { + if ( !k ) + { + errno = EDOM; + return 0.0; + } + } + + x = exp (t); + + if ( sign ) + { + px ^= 0x80000000; + SET_HIGH_WORD (x, px); + } + } + } + + return x; +} + +#endif _DOUBLE_IS_32BITS Index: Makefile.am =================================================================== --- Makefile.am (nonexistent) +++ Makefile.am (revision 1765) @@ -0,0 +1,195 @@ +## Process this file with automake to generate Makefile.in + +AUTOMAKE_OPTIONS = cygnus + +INCLUDES = -I$(srcdir)/../common $(NEWLIB_CFLAGS) $(CROSS_CFLAGS) $(TARGET_CFLAGS) + +src = s_acos.c s_frexp.c s_mathcnst.c \ + s_cos.c s_sinh.c \ + s_asin.c\ + s_asine.c s_cosh.c s_ispos.c s_numtest.c s_sqrt.c \ + s_exp.c s_ldexp.c s_pow.c s_tan.c \ + s_atan.c \ + s_atan2.c s_fabs.c s_log.c s_tanh.c \ + s_log10.c s_sin.c \ + s_floor.c s_sine.c \ + s_atangent.c s_logarithm.c \ + s_sineh.c \ + s_ceil.c s_isnan.c s_isinf.c \ + e_acosh.c e_atanh.c e_remainder.c \ + er_gamma.c er_lgamma.c \ + s_erf.c e_j0.c e_j1.c w_jn.c e_hypot.c \ + w_cabs.c w_drem.c s_asinh.c s_fmod.c \ + e_scalb.c s_infconst.c s_signif.c + +fsrc = sf_ceil.c \ + sf_acos.c sf_frexp.c \ + sf_cos.c sf_sinh.c \ + sf_asine.c sf_cosh.c sf_ispos.c sf_numtest.c sf_sqrt.c \ + sf_asin.c \ + sf_exp.c sf_ldexp.c sf_pow.c sf_tan.c \ + sf_atan2.c sf_fabs.c sf_tanh.c \ + sf_atan.c sf_log10.c sf_sin.c\ + sf_floor.c sf_sine.c \ + sf_atangent.c sf_logarithm.c sf_sineh.c \ + sf_log.c sf_sineh.c \ + sf_isnan.c sf_isinf.c \ + ef_acosh.c ef_atanh.c ef_remainder.c \ + erf_gamma.c erf_lgamma.c \ + sf_erf.c ef_j0.c ef_j1.c wf_jn.c ef_hypot.c \ + wf_cabs.c wf_drem.c sf_asinh.c sf_fmod.c \ + ef_scalb.c sf_signif.c + +libmathfp_la_LDFLAGS = -Xcompiler -nostdlib + +if USE_LIBTOOL +noinst_LTLIBRARIES = libmathfp.la +libmathfp_la_SOURCES = $(src) $(fsrc) +noinst_DATA = objectlist.awk.in +else +noinst_LIBRARIES = lib.a +lib_a_SOURCES = $(src) $(fsrc) +noinst_DATA = +endif # USE_LIBTOOL + +include $(srcdir)/../../Makefile.shared + +chobj = eacosh.def \ + eatanh.def \ + ehypot.def \ + eremainder.def \ + erlgamma.def \ + sacos.def \ + sasine.def \ + sasinh.def \ + satan.def \ + satan2.def \ + satangent.def \ + scosh.def \ + serf.def \ + sexp.def \ + sfabs.def \ + sfloor.def \ + sfmod.def \ + sfrexp.def \ + sisnan.def \ + sldexp.def \ + slog10.def \ + slogarithm.def \ + spow.def \ + ssine.def \ + ssineh.def \ + ssqrt.def \ + stan.def \ + stanh.def \ + wjn.def + +SUFFIXES = .def + +CHEW = ../../doc/makedoc -f $(srcdir)/../../doc/doc.str + +.c.def: + $(CHEW) < $< > $*.def 2> $*.ref + touch stmp-def + +TARGETDOC = ../tmp.texi + +doc: $(chobj) + cat $(srcdir)/mathfp.tex >> $(TARGETDOC) + +CLEANFILES = $(chobj) *.ref + +# Texinfo does not appear to support underscores in file names, so we +# name the .def files without underscores. + +eacosh.def: e_acosh.c + $(CHEW) < $(srcdir)/e_acosh.c >$@ 2>/dev/null + touch stmp-def +eatanh.def: e_atanh.c + $(CHEW) < $(srcdir)/e_atanh.c >$@ 2>/dev/null + touch stmp-def +ehypot.def: e_hypot.c + $(CHEW) < $(srcdir)/e_hypot.c >$@ 2>/dev/null + touch stmp-def +eremainder.def: e_remainder.c + $(CHEW) < $(srcdir)/e_remainder.c >$@ 2>/dev/null + touch stmp-def +erlgamma.def: er_lgamma.c + $(CHEW) < $(srcdir)/er_lgamma.c >$@ 2>/dev/null + touch stmp-def +sacos.def: s_acos.c + $(CHEW) < $(srcdir)/s_acos.c >$@ 2>/dev/null + touch stmp-def +sasine.def: s_asine.c + $(CHEW) < $(srcdir)/s_asine.c >$@ 2>/dev/null + touch stmp-def +sasinh.def: s_asinh.c + $(CHEW) < $(srcdir)/s_asinh.c >$@ 2>/dev/null + touch stmp-def +satan.def: s_atan.c + $(CHEW) < $(srcdir)/s_atan.c >$@ 2>/dev/null + touch stmp-def +satan2.def: s_atan2.c + $(CHEW) < $(srcdir)/s_atan2.c >$@ 2>/dev/null + touch stmp-def +satangent.def: s_atangent.c + $(CHEW) < $(srcdir)/s_atangent.c >$@ 2>/dev/null + touch stmp-def +scosh.def: s_cosh.c + $(CHEW) < $(srcdir)/s_cosh.c >$@ 2>/dev/null + touch stmp-def +serf.def: s_erf.c + $(CHEW) < $(srcdir)/s_erf.c >$@ 2>/dev/null + touch stmp-def +sexp.def: s_exp.c + $(CHEW) < $(srcdir)/s_exp.c >$@ 2>/dev/null + touch stmp-def +sfabs.def: s_fabs.c + $(CHEW) < $(srcdir)/s_fabs.c >$@ 2>/dev/null + touch stmp-def +sfloor.def: s_floor.c + $(CHEW) < $(srcdir)/s_floor.c >$@ 2>/dev/null + touch stmp-def +sfmod.def: s_fmod.c + $(CHEW) < $(srcdir)/s_fmod.c >$@ 2>/dev/null + touch stmp-def +sfrexp.def: s_frexp.c + $(CHEW) < $(srcdir)/s_frexp.c >$@ 2>/dev/null + touch stmp-def +sisnan.def: s_isnan.c + $(CHEW) < $(srcdir)/s_isnan.c >$@ 2>/dev/null + touch stmp-def +sldexp.def: s_ldexp.c + $(CHEW) < $(srcdir)/s_ldexp.c >$@ 2>/dev/null + touch stmp-def +slog10.def: s_log10.c + $(CHEW) < $(srcdir)/s_log10.c >$@ 2>/dev/null + touch stmp-def +slogarithm.def: s_logarithm.c + $(CHEW) < $(srcdir)/s_logarithm.c >$@ 2>/dev/null + touch stmp-def +spow.def: s_pow.c + $(CHEW) < $(srcdir)/s_pow.c >$@ 2>/dev/null + touch stmp-def +ssine.def: s_sine.c + $(CHEW) < $(srcdir)/s_sine.c >$@ 2>/dev/null + touch stmp-def +ssineh.def: s_sineh.c + $(CHEW) < $(srcdir)/s_sineh.c >$@ 2>/dev/null + touch stmp-def +ssqrt.def: s_sqrt.c + $(CHEW) < $(srcdir)/s_sqrt.c >$@ 2>/dev/null + touch stmp-def +stan.def: s_tan.c + $(CHEW) < $(srcdir)/s_tan.c >$@ 2>/dev/null + touch stmp-def +stanh.def: s_tanh.c + $(CHEW) < $(srcdir)/s_tanh.c >$@ 2>/dev/null + touch stmp-def +wjn.def: w_jn.c + $(CHEW) < $(srcdir)/w_jn.c >$@ 2>/dev/null + touch stmp-def + +# A partial dependency list. + +$(lib_a_OBJECTS): $(srcdir)/../../libc/include/math.h $(srcdir)/../common/fdlibm.h Index: sf_atangent.c =================================================================== --- sf_atangent.c (nonexistent) +++ sf_atangent.c (revision 1765) @@ -0,0 +1,140 @@ + +/* @(#)z_atangentf.c 1.0 98/08/13 */ +/****************************************************************** + * The following routines are coded directly from the algorithms + * and coefficients given in "Software Manual for the Elementary + * Functions" by William J. Cody, Jr. and William Waite, Prentice + * Hall, 1980. + ******************************************************************/ +/****************************************************************** + * Arctangent + * + * Input: + * x - floating point value + * + * Output: + * arctangent of x + * + * Description: + * This routine calculates arctangents. + * + *****************************************************************/ + +#include +#include "fdlibm.h" +#include "zmath.h" + +static const float ROOT3 = 1.732050807; +static const float a[] = { 0.0, 0.523598775, 1.570796326, + 1.047197551 }; +static const float q[] = { 0.1412500740e+1 }; +static const float p[] = { -0.4708325141, -0.5090958253e-1 }; + +float +_DEFUN (atangentf, (float, float, float, int), + float x _AND + float v _AND + float u _AND + int arctan2) +{ + float f, g, R, P, Q, A, res; + int N; + int branch = 0; + int expv, expu; + + /* Preparation for calculating arctan2. */ + if (arctan2) + { + if (u == 0.0) + if (v == 0.0) + { + errno = ERANGE; + return (z_notanum_f.f); + } + else + { + branch = 1; + res = __PI_OVER_TWO; + } + + if (!branch) + { + int e; + /* Get the exponent values of the inputs. */ + g = frexpf (v, &expv); + g = frexpf (u, &expu); + + /* See if a divide will overflow. */ + e = expv - expu; + if (e > FLT_MAX_EXP) + { + branch = 1; + res = __PI_OVER_TWO; + } + + /* Also check for underflow. */ + else if (e < FLT_MIN_EXP) + { + branch = 2; + res = 0.0; + } + } + } + + if (!branch) + { + if (arctan2) + f = fabsf (v / u); + else + f = fabsf (x); + + if (f > 1.0) + { + f = 1.0 / f; + N = 2; + } + else + N = 0; + + if (f > (2.0 - ROOT3)) + { + A = ROOT3 - 1.0; + f = (((A * f - 0.5) - 0.5) + f) / (ROOT3 + f); + N++; + } + + /* Check for values that are too small. */ + if (-z_rooteps_f < f && f < z_rooteps_f) + res = f; + + /* Calculate the Taylor series. */ + else + { + g = f * f; + P = (p[1] * g + p[0]) * g; + Q = g + q[0]; + R = P / Q; + + res = f + f * R; + } + + if (N > 1) + res = -res; + + res += a[N]; + } + + if (arctan2) + { + if (u < 0.0 || branch == 2) + res = __PI - res; + if (v < 0.0 || branch == 1) + res = -res; + } + else if (x < 0.0) + { + res = -res; + } + + return (res); +} Index: sf_asinh.c =================================================================== --- sf_asinh.c (nonexistent) +++ sf_asinh.c (revision 1765) @@ -0,0 +1,66 @@ +/* sf_asinh.c -- float version of s_asinh.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "fdlibm.h" + +#ifdef __STDC__ +static const float +#else +static float +#endif +one = 1.0000000000e+00, /* 0x3F800000 */ +ln2 = 6.9314718246e-01, /* 0x3f317218 */ +huge= 1.0000000000e+30; + +#ifdef __STDC__ + float asinhf(float x) +#else + float asinhf(x) + float x; +#endif +{ + float t,w; + __int32_t hx,ix; + GET_FLOAT_WORD(hx,x); + ix = hx&0x7fffffff; + if(ix>=0x7f800000) return x+x; /* x is inf or NaN */ + if(ix< 0x31800000) { /* |x|<2**-28 */ + if(huge+x>one) return x; /* return x inexact except 0 */ + } + if(ix>0x4d800000) { /* |x| > 2**28 */ + w = logf(fabsf(x))+ln2; + } else if (ix>0x40000000) { /* 2**28 > |x| > 2.0 */ + t = fabsf(x); + w = logf((float)2.0*t+one/(sqrtf(x*x+one)+t)); + } else { /* 2.0 > |x| > 2**-28 */ + t = x*x; + w =log1pf(fabsf(x)+t/(one+sqrtf(one+t))); + } + if(hx>0) return w; else return -w; +} + +#ifdef _DOUBLE_IS_32BITS + +#ifdef __STDC__ + double asinh(double x) +#else + double asinh(x) + double x; +#endif +{ + return (double) asinhf((float) x); +} + +#endif /* defined(_DOUBLE_IS_32BITS) */ Index: erf_lgamma.c =================================================================== --- erf_lgamma.c (nonexistent) +++ erf_lgamma.c (revision 1765) @@ -0,0 +1,244 @@ +/* erf_lgamma.c -- float version of er_lgamma.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + */ + +#include "fdlibm.h" + +#ifdef __STDC__ +static const float +#else +static float +#endif +two23= 8.3886080000e+06, /* 0x4b000000 */ +half= 5.0000000000e-01, /* 0x3f000000 */ +one = 1.0000000000e+00, /* 0x3f800000 */ +pi = 3.1415927410e+00, /* 0x40490fdb */ +a0 = 7.7215664089e-02, /* 0x3d9e233f */ +a1 = 3.2246702909e-01, /* 0x3ea51a66 */ +a2 = 6.7352302372e-02, /* 0x3d89f001 */ +a3 = 2.0580807701e-02, /* 0x3ca89915 */ +a4 = 7.3855509982e-03, /* 0x3bf2027e */ +a5 = 2.8905137442e-03, /* 0x3b3d6ec6 */ +a6 = 1.1927076848e-03, /* 0x3a9c54a1 */ +a7 = 5.1006977446e-04, /* 0x3a05b634 */ +a8 = 2.2086278477e-04, /* 0x39679767 */ +a9 = 1.0801156895e-04, /* 0x38e28445 */ +a10 = 2.5214456400e-05, /* 0x37d383a2 */ +a11 = 4.4864096708e-05, /* 0x383c2c75 */ +tc = 1.4616321325e+00, /* 0x3fbb16c3 */ +tf = -1.2148628384e-01, /* 0xbdf8cdcd */ +/* tt = -(tail of tf) */ +tt = 6.6971006518e-09, /* 0x31e61c52 */ +t0 = 4.8383611441e-01, /* 0x3ef7b95e */ +t1 = -1.4758771658e-01, /* 0xbe17213c */ +t2 = 6.4624942839e-02, /* 0x3d845a15 */ +t3 = -3.2788541168e-02, /* 0xbd064d47 */ +t4 = 1.7970675603e-02, /* 0x3c93373d */ +t5 = -1.0314224288e-02, /* 0xbc28fcfe */ +t6 = 6.1005386524e-03, /* 0x3bc7e707 */ +t7 = -3.6845202558e-03, /* 0xbb7177fe */ +t8 = 2.2596477065e-03, /* 0x3b141699 */ +t9 = -1.4034647029e-03, /* 0xbab7f476 */ +t10 = 8.8108185446e-04, /* 0x3a66f867 */ +t11 = -5.3859531181e-04, /* 0xba0d3085 */ +t12 = 3.1563205994e-04, /* 0x39a57b6b */ +t13 = -3.1275415677e-04, /* 0xb9a3f927 */ +t14 = 3.3552918467e-04, /* 0x39afe9f7 */ +u0 = -7.7215664089e-02, /* 0xbd9e233f */ +u1 = 6.3282704353e-01, /* 0x3f2200f4 */ +u2 = 1.4549225569e+00, /* 0x3fba3ae7 */ +u3 = 9.7771751881e-01, /* 0x3f7a4bb2 */ +u4 = 2.2896373272e-01, /* 0x3e6a7578 */ +u5 = 1.3381091878e-02, /* 0x3c5b3c5e */ +v1 = 2.4559779167e+00, /* 0x401d2ebe */ +v2 = 2.1284897327e+00, /* 0x4008392d */ +v3 = 7.6928514242e-01, /* 0x3f44efdf */ +v4 = 1.0422264785e-01, /* 0x3dd572af */ +v5 = 3.2170924824e-03, /* 0x3b52d5db */ +s0 = -7.7215664089e-02, /* 0xbd9e233f */ +s1 = 2.1498242021e-01, /* 0x3e5c245a */ +s2 = 3.2577878237e-01, /* 0x3ea6cc7a */ +s3 = 1.4635047317e-01, /* 0x3e15dce6 */ +s4 = 2.6642270386e-02, /* 0x3cda40e4 */ +s5 = 1.8402845599e-03, /* 0x3af135b4 */ +s6 = 3.1947532989e-05, /* 0x3805ff67 */ +r1 = 1.3920053244e+00, /* 0x3fb22d3b */ +r2 = 7.2193557024e-01, /* 0x3f38d0c5 */ +r3 = 1.7193385959e-01, /* 0x3e300f6e */ +r4 = 1.8645919859e-02, /* 0x3c98bf54 */ +r5 = 7.7794247773e-04, /* 0x3a4beed6 */ +r6 = 7.3266842264e-06, /* 0x36f5d7bd */ +w0 = 4.1893854737e-01, /* 0x3ed67f1d */ +w1 = 8.3333335817e-02, /* 0x3daaaaab */ +w2 = -2.7777778450e-03, /* 0xbb360b61 */ +w3 = 7.9365057172e-04, /* 0x3a500cfd */ +w4 = -5.9518753551e-04, /* 0xba1c065c */ +w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */ +w6 = -1.6309292987e-03; /* 0xbad5c4e8 */ + +#ifdef __STDC__ +static const float zero= 0.0000000000e+00; +#else +static float zero= 0.0000000000e+00; +#endif + +#ifdef __STDC__ + static float sin_pif(float x) +#else + static float sin_pif(x) + float x; +#endif +{ + float y,z; + __int32_t n,ix; + + GET_FLOAT_WORD(ix,x); + ix &= 0x7fffffff; + + if(ix<0x3e800000) return __kernel_sinf(pi*x,zero,0); + y = -x; /* x is assume negative */ + + /* + * argument reduction, make sure inexact flag not raised if input + * is an integer + */ + z = floorf(y); + if(z!=y) { /* inexact anyway */ + y *= (float)0.5; + y = (float)2.0*(y - floorf(y)); /* y = |x| mod 2.0 */ + n = (__int32_t) (y*(float)4.0); + } else { + if(ix>=0x4b800000) { + y = zero; n = 0; /* y must be even */ + } else { + if(ix<0x4b000000) z = y+two23; /* exact */ + GET_FLOAT_WORD(n,z); + n &= 1; + y = n; + n<<= 2; + } + } + switch (n) { + case 0: y = __kernel_sinf(pi*y,zero,0); break; + case 1: + case 2: y = __kernel_cosf(pi*((float)0.5-y),zero); break; + case 3: + case 4: y = __kernel_sinf(pi*(one-y),zero,0); break; + case 5: + case 6: y = -__kernel_cosf(pi*(y-(float)1.5),zero); break; + default: y = __kernel_sinf(pi*(y-(float)2.0),zero,0); break; + } + return -y; +} + + +#ifdef __STDC__ + float lgammaf_r(float x, int *signgamp) +#else + float lgammaf_r(x,signgamp) + float x; int *signgamp; +#endif +{ + float t,y,z,nadj,p,p1,p2,p3,q,r,w; + __int32_t i,hx,ix; + + GET_FLOAT_WORD(hx,x); + + /* purge off +-inf, NaN, +-0, and negative arguments */ + *signgamp = 1; + ix = hx&0x7fffffff; + if(ix>=0x7f800000) return x*x; + if(ix==0) return one/zero; + if(ix<0x1c800000) { /* |x|<2**-70, return -log(|x|) */ + if(hx<0) { + *signgamp = -1; + return -logf(-x); + } else return -logf(x); + } + if(hx<0) { + if(ix>=0x4b000000) /* |x|>=2**23, must be -integer */ + return one/zero; + t = sin_pif(x); + if(t==zero) return one/zero; /* -integer */ + nadj = logf(pi/fabsf(t*x)); + if(t=0x3f3b4a20) {y = one-x; i= 0;} + else if(ix>=0x3e6d3308) {y= x-(tc-one); i=1;} + else {y = x; i=2;} + } else { + r = zero; + if(ix>=0x3fdda618) {y=(float)2.0-x;i=0;} /* [1.7316,2] */ + else if(ix>=0x3F9da620) {y=x-tc;i=1;} /* [1.23,1.73] */ + else {y=x-one;i=2;} + } + switch(i) { + case 0: + z = y*y; + p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10)))); + p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11))))); + p = y*p1+p2; + r += (p-(float)0.5*y); break; + case 1: + z = y*y; + w = z*y; + p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */ + p2 = t1+w*(t4+w*(t7+w*(t10+w*t13))); + p3 = t2+w*(t5+w*(t8+w*(t11+w*t14))); + p = z*p1-(tt-w*(p2+y*p3)); + r += (tf + p); break; + case 2: + p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5))))); + p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5)))); + r += (-(float)0.5*y + p1/p2); + } + } + else if(ix<0x41000000) { /* x < 8.0 */ + i = (__int32_t)x; + t = zero; + y = x-(float)i; + p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); + q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); + r = half*y+p/q; + z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ + switch(i) { + case 7: z *= (y+(float)6.0); /* FALLTHRU */ + case 6: z *= (y+(float)5.0); /* FALLTHRU */ + case 5: z *= (y+(float)4.0); /* FALLTHRU */ + case 4: z *= (y+(float)3.0); /* FALLTHRU */ + case 3: z *= (y+(float)2.0); /* FALLTHRU */ + r += logf(z); break; + } + /* 8.0 <= x < 2**58 */ + } else if (ix < 0x5c800000) { + t = logf(x); + z = one/x; + y = z*z; + w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); + r = (x-half)*(t-one)+w; + } else + /* 2**58 <= x <= inf */ + r = x*(logf(x)-one); + if(hx<0) r = nadj - r; + return r; +}

Browse

Tools

Subversion Repositories or1k

Compare Revisions

Rev 1010 → Rev 1765