From 11efd5076a5e300a96514294a7e39f0b295553e1 Mon Sep 17 00:00:00 2001 From: Oliver Lehmann Date: Thu, 13 Nov 1986 21:40:52 +0000 Subject: [PATCH] add LCC - a better C-Compiler for WEGA --- WEGA/contrib/lcc/lib/libm/asin.c | 90 ++++++++++++ WEGA/contrib/lcc/lib/libm/atan.c | 105 ++++++++++++++ WEGA/contrib/lcc/lib/libm/erf.c | 107 ++++++++++++++ WEGA/contrib/lcc/lib/libm/exp.c | 68 +++++++++ WEGA/contrib/lcc/lib/libm/fabs.c | 13 ++ WEGA/contrib/lcc/lib/libm/floor.c | 28 ++++ WEGA/contrib/lcc/lib/libm/fmod.c | 32 ++++ WEGA/contrib/lcc/lib/libm/gamma.c | 128 ++++++++++++++++ WEGA/contrib/lcc/lib/libm/hypot.c | 48 ++++++ WEGA/contrib/lcc/lib/libm/j0.c | 217 ++++++++++++++++++++++++++++ WEGA/contrib/lcc/lib/libm/j1.c | 217 ++++++++++++++++++++++++++++ WEGA/contrib/lcc/lib/libm/jn.c | 143 ++++++++++++++++++ WEGA/contrib/lcc/lib/libm/log.c | 88 +++++++++++ WEGA/contrib/lcc/lib/libm/matherr.c | 13 ++ WEGA/contrib/lcc/lib/libm/pow.c | 89 ++++++++++++ WEGA/contrib/lcc/lib/libm/sin.c | 106 ++++++++++++++ WEGA/contrib/lcc/lib/libm/sinh.c | 91 ++++++++++++ WEGA/contrib/lcc/lib/libm/sqrt.c | 48 ++++++ WEGA/contrib/lcc/lib/libm/tan.c | 90 ++++++++++++ WEGA/contrib/lcc/lib/libm/tanh.c | 47 ++++++ 20 files changed, 1768 insertions(+) create mode 100644 WEGA/contrib/lcc/lib/libm/asin.c create mode 100644 WEGA/contrib/lcc/lib/libm/atan.c create mode 100644 WEGA/contrib/lcc/lib/libm/erf.c create mode 100644 WEGA/contrib/lcc/lib/libm/exp.c create mode 100644 WEGA/contrib/lcc/lib/libm/fabs.c create mode 100644 WEGA/contrib/lcc/lib/libm/floor.c create mode 100644 WEGA/contrib/lcc/lib/libm/fmod.c create mode 100644 WEGA/contrib/lcc/lib/libm/gamma.c create mode 100644 WEGA/contrib/lcc/lib/libm/hypot.c create mode 100644 WEGA/contrib/lcc/lib/libm/j0.c create mode 100644 WEGA/contrib/lcc/lib/libm/j1.c create mode 100644 WEGA/contrib/lcc/lib/libm/jn.c create mode 100644 WEGA/contrib/lcc/lib/libm/log.c create mode 100644 WEGA/contrib/lcc/lib/libm/matherr.c create mode 100644 WEGA/contrib/lcc/lib/libm/pow.c create mode 100644 WEGA/contrib/lcc/lib/libm/sin.c create mode 100644 WEGA/contrib/lcc/lib/libm/sinh.c create mode 100644 WEGA/contrib/lcc/lib/libm/sqrt.c create mode 100644 WEGA/contrib/lcc/lib/libm/tan.c create mode 100644 WEGA/contrib/lcc/lib/libm/tanh.c diff --git a/WEGA/contrib/lcc/lib/libm/asin.c b/WEGA/contrib/lcc/lib/libm/asin.c new file mode 100644 index 0000000..76b7d78 --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/asin.c @@ -0,0 +1,90 @@ +/* ~(#)asin.c 2.1 */ +/* ~(#)asin.c 1.12 */ +/*LINTLIBRARY*/ +/* + * C program for double-precision asin/acos. + * Returns EDOM error and value 0 if |argument| > 1. + * Algorithm and coefficients from Cody and Waite (1980). + * Calls sqrt if |argument| > 0.5. + */ + +#include +#include +#include + +double +asin(x) +double x; +{ + extern double asin_acos(); + + return (asin_acos(x, 0)); +} + +static double +asin_acos(x, acosflag) +register double x; +int acosflag; +{ + register double y; + register int neg = 0, large = 0; + struct exception exc; + + exc.arg1 = x; + if (x < 0) { + x = -x; + neg++; + } + if (x > 1) { + exc.type = DOMAIN; + exc.name = acosflag ? "acos" : "asin"; + exc.retval = 0.0; + if (!matherr(&exc)) { + (void) write(2, exc.name, 4); + (void) write(2, ": DOMAIN error\n", 15); + errno = EDOM; + } + return (exc.retval); + } + if (x > X_EPS) { /* skip for efficiency and to prevent underflow */ + static double p[] = { + -0.69674573447350646411e0, + 0.10152522233806463645e2, + -0.39688862997504877339e2, + 0.57208227877891731407e2, + -0.27368494524164255994e2, + }, q[] = { + 1.0, + -0.23823859153670238830e2, + 0.15095270841030604719e3, + -0.38186303361750149284e3, + 0.41714430248260412556e3, + -0.16421096714498560795e3, + }; + + if (x <= 0.5) + y = x * x; + else { + large++; + y = 0.5 - 0.5 * x; + x = -sqrt(y); + x += x; + } + x += x * y * _POLY4(y, p)/_POLY5(y, q); + } + if (acosflag) { + if (!neg) + x = -x; + return (!large ? M_PI_2 + x : neg ? M_PI + x : x); + } + if (large) + x += M_PI_2; + return (neg ? -x : x); +} + +double +acos(x) +double x; +{ + return (asin_acos(x, 1)); +} diff --git a/WEGA/contrib/lcc/lib/libm/atan.c b/WEGA/contrib/lcc/lib/libm/atan.c new file mode 100644 index 0000000..8a33042 --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/atan.c @@ -0,0 +1,105 @@ +/* ~(#)atan.c 2.1 */ +/* ~(#)atan.c 1.10 */ +/*LINTLIBRARY*/ +/* + * atan returns the arctangent of its double-precision argument, + * in the range [-pi/2, pi/2]. + * There are no error returns. + * + * atan2(y, x) returns the arctangent of y/x, + * in the range [-pi, pi]. + * atan2 discovers what quadrant the angle is in and calls atan. + * atan2 returns EDOM error and value 0 if both arguments are zero. + * + * Coefficients are #5077 from Hart & Cheney (19.56D). + */ + +#include +#include +#include +#define SQ2_M_1 0.41421356237309504880 +#define SQ2_P_1 2.41421356237309504880 +#define NEG 1 +#define INV 2 +#define MID 4 + +double +atan(x) +register double x; +{ + register int type = 0; + + if (x < 0) { + x = -x; + type = NEG; + } + if (x > SQ2_P_1) { + x = 1/x; + type |= INV; + } else if (x > SQ2_M_1) { + x = (x - 1)/(x + 1); + type |= MID; + } + if (x < -X_EPS || x > X_EPS) { + /* skip for efficiency and to prevent underflow */ + static double p[] = { + 0.161536412982230228262e2, + 0.26842548195503973794141e3, + 0.11530293515404850115428136e4, + 0.178040631643319697105464587e4, + 0.89678597403663861959987488e3, + }, q[] = { + 1.0, + 0.5895697050844462222791e2, + 0.536265374031215315104235e3, + 0.16667838148816337184521798e4, + 0.207933497444540981287275926e4, + 0.89678597403663861962481162e3, + }; + register double xsq = x * x; + + x *= _POLY4(xsq, p)/_POLY5(xsq, q); + } + if (type & INV) + x = M_PI_2 - x; + else if (type & MID) + x += M_PI_4; + return (type & NEG ? -x : x); +} + +double +atan2(y, x) +register double y, x; +{ + register int neg_y = 0; + struct exception exc; + + exc.arg1 = y; + if (!x && !y) { + exc.type = DOMAIN; + exc.name = "atan2"; + exc.arg2 = x; + exc.retval = 0.0; + if (!matherr(&exc)) { + (void) write(2, "atan2: DOMAIN error\n", 20); + errno = EDOM; + } + return (exc.retval); + } + /* + * The next lines determine if |x| is negligible compared to |y|, + * without dividing, and without adding values of the same sign. + */ + if (exc.arg1 < 0) { + exc.arg1 = -exc.arg1; + neg_y++; + } + if (exc.arg1 - _ABS(x) == exc.arg1) + return (neg_y ? -M_PI_2 : M_PI_2); + /* + * The next line assumes that if y/x underflows the result + * is zero with no error indication, so it's safe to divide. + */ + return ((y = atan(y/x)) == 0 || x > 0 ? y : + neg_y ? y - M_PI : y + M_PI); +} diff --git a/WEGA/contrib/lcc/lib/libm/erf.c b/WEGA/contrib/lcc/lib/libm/erf.c new file mode 100644 index 0000000..6d1f987 --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/erf.c @@ -0,0 +1,107 @@ +/* ~(#)erf.c 2.1 */ +/* ~(#)erf.c 1.10 */ +/*LINTLIBRARY*/ +/* + * erf returns the error function of its double-precision argument. + * erfc(x) returns 1 - erf(x). + * + * erf(x) is defined by + * ${2 over sqrt pi} int from 0 to x e sup {- t sup 2} dt$. + * + * The entry for erfc is provided because of the + * extreme loss of relative accuracy if erf(x) is + * called for large x and the result subtracted + * from 1 (e.g. for x = 5, 12 places are lost). + * + * There are no error returns. + * + * Calls exp for |x| > 0.5. + * + * Coefficients for large x are #5667 from Hart & Cheney (18.72D). + */ + +#include +/* approx sqrt(log(MAXDOUBLE)) */ +#if u3b || u3b5 +#define MAXVAL 27.23 +#else +#define MAXVAL 9.27 +#endif +#define DPOLYD(y, p, q) for (n = d = 0, i = sizeof(p)/sizeof(p[0]); --i >= 0; ) \ + { n = n * y + p[i]; d = d * y + q[i]; } + +static double p1[] = { + 0.804373630960840172832162e5, + 0.740407142710151470082064e4, + 0.301782788536507577809226e4, + 0.380140318123903008244444e2, + 0.143383842191748205576712e2, + -.288805137207594084924010e0, + 0.007547728033418631287834e0, +}, q1[] = { + 0.804373630960840172826266e5, + 0.342165257924628539769006e5, + 0.637960017324428279487120e4, + 0.658070155459240506326937e3, + 0.380190713951939403753468e2, + 1.0, + 0.0, +}; +static double p2[] = { + 0.18263348842295112592168999e4, + 0.28980293292167655611275846e4, + 0.2320439590251635247384768711e4, + 0.1143262070703886173606073338e4, + 0.3685196154710010637133875746e3, + 0.7708161730368428609781633646e2, + 0.9675807882987265400604202961e1, + 0.5641877825507397413087057563e0, + 0.0, +}, q2[] = { + 0.18263348842295112595576438e4, + 0.495882756472114071495438422e4, + 0.60895424232724435504633068e4, + 0.4429612803883682726711528526e4, + 0.2094384367789539593790281779e4, + 0.6617361207107653469211984771e3, + 0.1371255960500622202878443578e3, + 0.1714980943627607849376131193e2, + 1.0, +}; + +double +erf(x) +register double x; +{ + int neg = 0; + + if (x < 0) { + x = -x; + neg++; + } + if (x > 0.5) + x = 1 - erfc(x); + else { + register double n, d, xsq = x * x; + register int i; + + DPOLYD(xsq, p1, q1); + x *= M_2_SQRTPI * n/d; + } + return (neg ? -x : x); +} + +double +erfc(x) +register double x; +{ + register double n, d; + register int i; + + if (x < 0.5) + return (1 - erf(x)); + if (x >= MAXVAL) /* exp(-x * x) sure to underflow */ + return (0.0); + DPOLYD(x, p2, q2); + return (exp(-x * x) * n/d); +} diff --git a/WEGA/contrib/lcc/lib/libm/exp.c b/WEGA/contrib/lcc/lib/libm/exp.c new file mode 100644 index 0000000..d4ae10f --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/exp.c @@ -0,0 +1,68 @@ +/* ~(#)exp.c 2.1 */ +/* ~(#)exp.c 1.11 */ +/*LINTLIBRARY*/ +/* + * exp returns the exponential function of its double-precision argument. + * Returns ERANGE error and value 0 if argument too small, + * value HUGE if argument too large. + * Algorithm and coefficients from Cody and Waite (1980). + * Calls ldexp. + */ + +#include +#include +#include + +double +exp(x) +register double x; +{ + static double p[] = { + 0.31555192765684646356e-4, + 0.75753180159422776666e-2, + 0.25000000000000000000e0, + }, q[] = { + 0.75104028399870046114e-6, + 0.63121894374398503557e-3, + 0.56817302698551221787e-1, + 0.50000000000000000000e0, + }; + register double y; + register int n; + struct exception exc; + + exc.arg1 = x; + if (x < 0) + x = -x; + if (x < X_EPS) /* use first term of Taylor series expansion */ + return (1.0 + exc.arg1); + exc.name = "exp"; + if (exc.arg1 <= LN_MINDOUBLE) { + if (exc.arg1 == LN_MINDOUBLE) /* protect against roundoff */ + return (MINDOUBLE); /* causing ldexp to underflow */ + exc.type = UNDERFLOW; + exc.retval = 0.0; + if (!matherr(&exc)) + errno = ERANGE; + return (exc.retval); + } + if (exc.arg1 >= LN_MAXDOUBLE) { + if (exc.arg1 == LN_MAXDOUBLE) /* protect against roundoff */ + return (MAXDOUBLE); /* causing ldexp to overflow */ + exc.type = OVERFLOW; + exc.retval = HUGE; + if (!matherr(&exc)) + errno = ERANGE; + return (exc.retval); + } + n = (int)(x * M_LOG2E + 0.5); + y = (double)n; + _REDUCE(int, x, y, 0.693359375, -2.1219444005469058277e-4); + if (exc.arg1 < 0) { + x = -x; + n = -n; + } + y = x * x; + x *= _POLY2(y, p); + return (ldexp(0.5 + x/(_POLY3(y, q) - x), n + 1)); +} diff --git a/WEGA/contrib/lcc/lib/libm/fabs.c b/WEGA/contrib/lcc/lib/libm/fabs.c new file mode 100644 index 0000000..d1bb12c --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/fabs.c @@ -0,0 +1,13 @@ +/* ~(#)fabs.c 2.1 */ +/* ~(#)fabs.c 1.4 */ +/*LINTLIBRARY*/ +/* + * fabs returns the absolute value of its double-precision argument. + */ + +double +fabs(x) +register double x; +{ + return (x < 0 ? -x : x); +} diff --git a/WEGA/contrib/lcc/lib/libm/floor.c b/WEGA/contrib/lcc/lib/libm/floor.c new file mode 100644 index 0000000..762a6ea --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/floor.c @@ -0,0 +1,28 @@ +/* ~(#)floor.c 2.1 */ +/* ~(#)floor.c 1.6 */ +/*LINTLIBRARY*/ +/* + * floor(x) returns the largest integer (as a double-precision number) + * not greater than x. + * ceil(x) returns the smallest integer not less than x. + */ + +extern double modf(); + +double +floor(x) +double x; +{ + double y; /* can't be in register because of modf() below */ + + return (modf(x, &y) < 0 ? y - 1 : y); +} + +double +ceil(x) +double x; +{ + double y; /* can't be in register because of modf() below */ + + return (modf(x, &y) > 0 ? y + 1 : y); +} diff --git a/WEGA/contrib/lcc/lib/libm/fmod.c b/WEGA/contrib/lcc/lib/libm/fmod.c new file mode 100644 index 0000000..4e4e672 --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/fmod.c @@ -0,0 +1,32 @@ +/* ~(#)fmod.c 2.1.1.1 */ +/* ~(#)fmod.c 1.8 */ +/*LINTLIBRARY*/ +/* + * fmod(x, y) returns the remainder of x on division by y, + * with the same sign as x, + * except that if |y| << |x|, it returns 0. + */ + +#include + +double +fmod(x, y) +register double x, y; +{ + double d; /* can't be in register because of modf() below */ + + /* + * The next lines determine if |y| is negligible compared to |x|, + * without dividing, and without adding values of the same sign. + */ + if( y==0.0 ) return( x ); /* fmod(x,0.0)==x */ + d = _ABS(x); + if (d - _ABS(y) == d) + return (0.0); +#if !defined(pdp11) && !defined(u3b5) /* pdp11 and u3b5 "cc" cant handle + cast of a double to a void */ + (void) +#endif + modf(x/y, &d); /* now it's safe to divide without overflow */ + return (x - d * y); +} diff --git a/WEGA/contrib/lcc/lib/libm/gamma.c b/WEGA/contrib/lcc/lib/libm/gamma.c new file mode 100644 index 0000000..38c5113 --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/gamma.c @@ -0,0 +1,128 @@ +/* ~(#)gamma.c 2.1.1.2 */ +/* ~(#)gamma.c 1.14 */ +/*LINTLIBRARY*/ +/* + * gamma returns the log of the absolute value of the gamma function + * of its double-precision argument. + * The sign of the gamma function is returned in the + * external integer signgam. + * Returns EDOM error and value HUGE if argument is non-negative integer. + * Returns ERANGE error and value HUGE if the correct value would overflow. + * + * The coefficients for expansion around zero + * are #5243 from Hart & Cheney; for expansion + * around infinity they are #5404. + * + * Calls log and sin. + */ + +#include +#include +#include +#define X_MAX (3.0 * H_PREC) +#define GOOBIE 0.9189385332046727417803297 + +int signgam; + +double +gamma(x) +register double x; +{ + extern double pos_gamma(); + struct exception exc; + + exc.type = 0; + exc.name = "gamma"; + exc.arg1 = x; + exc.retval = HUGE; + signgam = 1; + if (x > 0) + x = pos_gamma(x, &exc); + else { + static double pi = M_PI; + double temp; /* can't be in register because of modf() below */ + + if (!modf(x = -x, &temp)) {/* DOMAIN if x is negative integer */ + exc.type = DOMAIN; + if (!matherr(&exc)) { + (void) write(2, "gamma: DOMAIN error\n", 20); + errno = EDOM; + } + return (exc.retval); + } + if (x >= X_MAX) + exc.type = OVERFLOW; + else { + if ((temp = sin(pi * x)) < 0) + temp = -temp; + else + signgam = -1; + return (-(log(x * temp/pi) + pos_gamma(x, &exc))); + } + } + if (exc.type != OVERFLOW) + return (x); + if (!matherr(&exc)) + errno = ERANGE; + return (exc.retval); +} + +static double +pos_gamma(x, excp) +register double x; +struct exception *excp; +{ + static double p2[] = { + -0.67449507245925289918e1, + -0.50108693752970953015e2, + -0.43933044406002567613e3, + -0.20085274013072791214e4, + -0.87627102978521489560e4, + -0.20886861789269887364e5, + -0.42353689509744089647e5, + }, q2[] = { + 1.0, + -0.23081551524580124562e2, + 0.18949823415702801641e3, + -0.49902852662143904834e3, + -0.15286072737795220248e4, + 0.99403074150827709015e4, + -0.29803853309256649932e4, + -0.42353689509744090010e5, + }; + register double y, z; + + if (x > 8) { /* asymptotic approximation */ + static double p[] = { + -0.1633436431e-2, + 0.83645878922e-3, + -0.5951896861197e-3, + 0.793650576493454e-3, + -0.277777777735865004e-2, + 0.83333333333333101837e-1, + }; + + if (x >= MAXDOUBLE/LN_MAXDOUBLE) { + excp->type = OVERFLOW; + return (excp->retval); + } + z = (x - 0.5) * log(x) - x + GOOBIE; + if (x > X_MAX) + return (z); + x = 1/x; + y = x * x; + return (z + x * _POLY5(y, p)); + } + y = 1; + if (x < y) + y /= (x * (y + x)); + else if (x < 2) { + y /= x; + x -= 1; + } else { + for ( ; x >= 3; y *= x) + x -= 1; + x -= 2; + } + return (log(y * _POLY6(x, p2)/_POLY7(x, q2))); +} diff --git a/WEGA/contrib/lcc/lib/libm/hypot.c b/WEGA/contrib/lcc/lib/libm/hypot.c new file mode 100644 index 0000000..91b083c --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/hypot.c @@ -0,0 +1,48 @@ +/* ~(#)hypot.c 2.1 */ +/* ~(#)hypot.c 1.10 */ +/*LINTLIBRARY*/ +/* + * hypot(a, b) returns sqrt(a^2 + b^2), avoiding unnecessary overflows. + * Returns ERANGE error and value HUGE if the correct value would overflow. + */ + +#include +#include +#include +#define ITERATIONS 4 + +double +hypot(a, b) +register double a, b; +{ + register double t; + register int i = ITERATIONS; + struct exception exc; + + if ((exc.arg1 = a) < 0) + a = -a; + if ((exc.arg2 = b) < 0) + b = -b; + if (a > b) { /* make sure |a| <= |b| */ + t = a; + a = b; + b = t; + } + if ((t = a/b) < X_EPS) /* t <= 1 */ + return (b); /* t << 1 */ + a = 1 + t * t; /* a = 1 + (a/b)^2 */ + t = 0.5 + 0.5 * a; /* first guess for sqrt */ + do { + t = 0.5 * (t + a/t); + } while (--i > 0); /* t <= sqrt(2) */ + if (b < MAXDOUBLE/M_SQRT2) /* result can't overflow */ + return (t * b); + if ((t *= 0.5 * b) < MAXDOUBLE/2) /* result won't overflow */ + return (t + t); + exc.type = OVERFLOW; + exc.name = "hypot"; + exc.retval = HUGE; + if (!matherr(&exc)) + errno = ERANGE; + return (exc.retval); +} diff --git a/WEGA/contrib/lcc/lib/libm/j0.c b/WEGA/contrib/lcc/lib/libm/j0.c new file mode 100644 index 0000000..2f57a84 --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/j0.c @@ -0,0 +1,217 @@ +/* ~(#)j0.c 2.1 */ +/* ~(#)j0.c 1.13 */ +/*LINTLIBRARY*/ +/* + * Double-precision Bessel's function + * of the first and second kinds + * of order zero. + * + * j0(x) returns the value of J0(x) + * for all real values of x. + * + * Returns ERANGE error and value 0 for large arguments. + * Calls sin, cos, sqrt. + * + * There is a niggling bug in J0 that + * causes errors up to 2e-16 for x in the + * interval [-8, 8]. + * The bug is caused by an inappropriate order + * of summation of the series. + * + * Coefficients are from Hart & Cheney. + * #5849 (19.22D) + * #6549 (19.25D) + * #6949 (19.41D) + * + * y0(x) returns the value of Y0(x) + * for positive real values of x. + * Returns EDOM error and value -HUGE if argument <= 0. + * + * Calls sin, cos, sqrt, log, j0. + * + * The values of Y0 have not been checked + * to more than ten places. + * + * Coefficients are from Hart & Cheney. + * #6245 (18.78D) + * #6549 (19.25D) + * #6949 (19.41D) + */ + +#include +#include +#include +#define P2_0_Q2_0 0.999999999999999999944688442 +#define P3_0_Q3_0 -0.0156249999999999999611615235 +#define P4_0_Q4_0 0.073804295108687225110222 +#define DPOLYD(y, p, q) for (n = d = 0, i = sizeof(p)/sizeof(p[0]); --i >= 0; ) \ + { n = n * y + p[i]; d = d * y + q[i]; } + +static double tpi = 0.6366197723675813430755350535; +static double p1[] = { + 0.4933787251794133561816813446e21, + -.1179157629107610536038440800e21, + 0.6382059341072356562289432465e19, + -.1367620353088171386865416609e18, + 0.1434354939140344111664316553e16, + -.8085222034853793871199468171e13, + 0.2507158285536881945555156435e11, + -.4050412371833132706360663322e8, + 0.2685786856980014981415848441e5, +}, q1[] = { + 0.4933787251794133562113278438e21, + 0.5428918384092285160200195092e19, + 0.3024635616709462698627330784e17, + 0.1127756739679798507056031594e15, + 0.3123043114941213172572469442e12, + 0.6699987672982239671814028660e9, + 0.1114636098462985378182402543e7, + 0.1363063652328970604442810507e4, + 1.0, +}; +static double p2[] = { + 0.5393485083869438325262122897e7, + 0.1233238476817638145232406055e8, + 0.8413041456550439208464315611e7, + 0.2016135283049983642487182349e7, + 0.1539826532623911470917825993e6, + 0.2485271928957404011288128951e4, + 0.0, +}, q2[] = { + 0.5393485083869438325560444960e7, + 0.1233831022786324960844856182e8, + 0.8426449050629797331554404810e7, + 0.2025066801570134013891035236e7, + 0.1560017276940030940592769933e6, + 0.2615700736920839685159081813e4, + 1.0, +}; +static double p3[] = { + -.3984617357595222463506790588e4, + -.1038141698748464093880530341e5, + -.8239066313485606568803548860e4, + -.2365956170779108192723612816e4, + -.2262630641933704113967255053e3, + -.4887199395841261531199129300e1, + 0.0, +}, q3[] = { + 0.2550155108860942382983170882e6, + 0.6667454239319826986004038103e6, + 0.5332913634216897168722255057e6, + 0.1560213206679291652539287109e6, + 0.1570489191515395519392882766e5, + 0.4087714673983499223402830260e3, + 1.0, +}; +static double p4[] = { + -.2750286678629109583701933175e20, + 0.6587473275719554925999402049e20, + -.5247065581112764941297350814e19, + 0.1375624316399344078571335453e18, + -.1648605817185729473122082537e16, + 0.1025520859686394284509167421e14, + -.3436371222979040378171030138e11, + 0.5915213465686889654273830069e8, + -.4137035497933148554125235152e5, +}, q4[] = { + 0.3726458838986165881989980e21, + 0.4192417043410839973904769661e19, + 0.2392883043499781857439356652e17, + 0.9162038034075185262489147968e14, + 0.2613065755041081249568482092e12, + 0.5795122640700729537480087915e9, + 0.1001702641288906265666651753e7, + 0.1282452772478993804176329391e4, + 1.0, +}; + +extern double j0_asympt(); + +double +j0(x) +register double x; +{ + register double n, d; + register int i; + + if ((n = x) < 0) + x = -x; + if (x > 8) + return (j0_asympt(x, n, 1)); + if (x < X_EPS) + return (1); + x *= x; + DPOLYD(x, p1, q1); + return (n/d); +} + +double +y0(x) +register double x; +{ + register double n, d, y, z; + register int i; + + if (x <= 0) { + struct exception exc; + + exc.type = DOMAIN; + exc.name = "y0"; + exc.arg1 = x; + exc.retval = -HUGE; + if (!matherr(&exc)) { + (void) write(2, "y0: DOMAIN error\n", 17); + errno = EDOM; + } + return (exc.retval); + } + if (x > 8) + return (j0_asympt(x, x, 0)); + y = tpi * log(x); + if (x < X_EPS) + return (y - P4_0_Q4_0); + z = x * x; + DPOLYD(z, p4, q4); + return (n/d + y * j0(x)); +} + +static double +j0_asympt(x, n, j0flag) +register double x, n; +int j0flag; +{ + register double z, d, pzero, qzero; + register int i; + + if (x > X_TLOSS) { + struct exception exc; + + exc.type = TLOSS; + exc.name = j0flag ? "j0" : "y0"; + exc.arg1 = n; + exc.retval = 0.0; + if (!matherr(&exc)) { + (void) write(2, exc.name, 2); + (void) write(2, ": TLOSS error\n", 14); + errno = ERANGE; + } + return (exc.retval); + } + if (x > X_PLOSS) { + pzero = P2_0_Q2_0; + qzero = P3_0_Q3_0; + } else { + z = 64/(x * x); + DPOLYD(z, p2, q2); + pzero = n/d; + DPOLYD(z, p3, q3); + qzero = n/d; + } + qzero *= 8/x; + z = sqrt(tpi/x); + pzero *= z; + qzero *= z; + x -= M_PI_4; + return (j0flag ? pzero * cos(x) - qzero * sin(x) : + pzero * sin(x) + qzero * cos(x)); +} diff --git a/WEGA/contrib/lcc/lib/libm/j1.c b/WEGA/contrib/lcc/lib/libm/j1.c new file mode 100644 index 0000000..569040b --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/j1.c @@ -0,0 +1,217 @@ +/* ~(#)j1.c 2.1 */ +/* ~(#)j1.c 1.11 */ +/*LINTLIBRARY*/ +/* + * Double-precision Bessel's function + * of the first and second kinds + * of order one. + * + * j1(x) returns the value of J1(x) + * for all real values of x. + * + * Returns ERANGE error and value 0 for large arguments. + * Calls sin, cos, sqrt. + * + * There is a niggling bug in J1 that + * causes errors up to 2e-16 for x in the + * interval [-8, 8]. + * The bug is caused by an inappropriate order + * of summation of the series. + * + * Coefficients are from Hart & Cheney. + * #6050 (20.98D) + * #6750 (19.19D) + * #7150 (19.35D) + * + * y1(x) returns the value of Y1(x) + * for positive real values of x. + * Returns EDOM error and value -HUGE if argument <= 0. + * + * Calls sin, cos, sqrt, log, j1. + * + * The values of Y1 have not been checked + * to more than ten places. + * + * Coefficients are from Hart & Cheney. + * #6447 (22.18D) + * #6750 (19.19D) + * #7150 (19.35D) + */ + +#include +#include +#include +#define P1_0_Q1_0 0.4999999999999999999989557017 +#define P2_0_Q2_0 1.0000000000000000000646346901 +#define P3_0_Q3_0 0.046874999999999999955398015174 +#define DPOLYD(y, p, q) for (n = d = 0, i = sizeof(p)/sizeof(p[0]); --i >= 0; ) \ + { n = n * y + p[i]; d = d * y + q[i]; } + +static double tpi = 0.6366197723675813430755350535; +static double p1[] = { + 0.581199354001606143928050809e21, + -.6672106568924916298020941484e20, + 0.2316433580634002297931815435e19, + -.3588817569910106050743641413e17, + 0.2908795263834775409737601689e15, + -.1322983480332126453125473247e13, + 0.3413234182301700539091292655e10, + -.4695753530642995859767162166e7, + 0.2701122710892323414856790990e4, +}, q1[] = { + 0.1162398708003212287858529400e22, + 0.1185770712190320999837113348e20, + 0.6092061398917521746105196863e17, + 0.2081661221307607351240184229e15, + 0.5243710262167649715406728642e12, + 0.1013863514358673989967045588e10, + 0.1501793594998585505921097578e7, + 0.1606931573481487801970916749e4, + 1.0, +}; +static double p2[] = { + -.4435757816794127857114720794e7, + -.9942246505077641195658377899e7, + -.6603373248364939109255245434e7, + -.1523529351181137383255105722e7, + -.1098240554345934672737413139e6, + -.1611616644324610116477412898e4, + 0.0, +}, q2[] = { + -.4435757816794127856828016962e7, + -.9934124389934585658967556309e7, + -.6585339479723087072826915069e7, + -.1511809506634160881644546358e7, + -.1072638599110382011903063867e6, + -.1455009440190496182453565068e4, + 1.0, +}; +static double p3[] = { + 0.3322091340985722351859704442e5, + 0.8514516067533570196555001171e5, + 0.6617883658127083517939992166e5, + 0.1849426287322386679652009819e5, + 0.1706375429020768002061283546e4, + 0.3526513384663603218592175580e2, + 0.0, +}, q3[] = { + 0.7087128194102874357377502472e6, + 0.1819458042243997298924553839e7, + 0.1419460669603720892855755253e7, + 0.4002944358226697511708610813e6, + 0.3789022974577220264142952256e5, + 0.8638367769604990967475517183e3, + 1.0, +}; +static double p4[] = { + -.9963753424306922225996744354e23, + 0.2655473831434854326894248968e23, + -.1212297555414509577913561535e22, + 0.2193107339917797592111427556e20, + -.1965887462722140658820322248e18, + 0.9569930239921683481121552788e15, + -.2580681702194450950541426399e13, + 0.3639488548124002058278999428e10, + -.2108847540133123652824139923e7, + 0.0, +}, q4[] = { + 0.5082067366941243245314424152e24, + 0.5435310377188854170800653097e22, + 0.2954987935897148674290758119e20, + 0.1082258259408819552553850180e18, + 0.2976632125647276729292742282e15, + 0.6465340881265275571961681500e12, + 0.1128686837169442121732366891e10, + 0.1563282754899580604737366452e7, + 0.1612361029677000859332072312e4, + 1.0, +}; + +extern double j1_asympt(); + +double +j1(x) +register double x; +{ + register double n, d, y; + register int i; + + if ((y = x) < 0) + x = -x; + if (x > 8) + return (j1_asympt(x, y, 1)); + if (x < X_EPS) + return (P1_0_Q1_0 * y); + x *= x; + DPOLYD(x, p1, q1); + return (y * n/d); +} + +double +y1(x) +register double x; +{ + register double n, d, z; + register int i; + + if (x <= 0) { + struct exception exc; + + exc.type = DOMAIN; + exc.name = "y1"; + exc.arg1 = x; + exc.retval = -HUGE; + if (!matherr(&exc)) { + (void) write(2, "y1: DOMAIN error\n", 17); + errno = EDOM; + } + return (exc.retval); + } + if (x > 8) + return (j1_asympt(x, x, 0)); + z = x * x; + DPOLYD(z, p4, q4); + return (x * n/d + tpi * (j1(x) * log(x) - 1/x)); +} + +static double +j1_asympt(x, n, j1flag) +register double x, n; +int j1flag; +{ + register double z, d, pzero, qzero; + register int i; + struct exception exc; + + exc.arg1 = n; + if (x > X_TLOSS) { + exc.type = TLOSS; + exc.name = j1flag ? "j1" : "y1"; + exc.retval = 0.0; + if (!matherr(&exc)) { + (void) write(2, exc.name, 2); + (void) write(2, ": TLOSS error\n", 14); + errno = ERANGE; + } + return (exc.retval); + } + if (x > X_PLOSS) { + pzero = P2_0_Q2_0; + qzero = P3_0_Q3_0; + } else { + z = 64/(x * x); + DPOLYD(z, p2, q2); + pzero = n/d; + DPOLYD(z, p3, q3); + qzero = n/d; + } + qzero *= 8/x; + z = sqrt(tpi/x); + pzero *= z; + qzero *= z; + x -= 3 * M_PI_4; + if (!j1flag) + return (pzero * sin(x) + qzero * cos(x)); + x = pzero * cos(x) - qzero * sin(x); + return (exc.arg1 < 0 ? -x : x); +} diff --git a/WEGA/contrib/lcc/lib/libm/jn.c b/WEGA/contrib/lcc/lib/libm/jn.c new file mode 100644 index 0000000..8092ead --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/jn.c @@ -0,0 +1,143 @@ +/* ~(#)jn.c 2.1 */ +/* ~(#)jn.c 1.10 */ +/*LINTLIBRARY*/ +/* + * Double-precision Bessel's function of + * the first and second kinds and of + * integer order. + * + * jn(n, x) returns the value of Jn(x) for all + * integer values of n and all real values + * of x. + * Returns ERANGE error and value 0 for large arguments. + * Calls j0, j1. + * + * For n = 0, j0(x) is called, + * For n = 1, j1(x) is called. + * For n < x, forward recursion is used starting + * from values of j0(x) and j1(x). + * For n > x, 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 y0 and y1 are called, and that forward recursion + * is used for values of n > 1. + * Returns EDOM error and value -HUGE if argument <= 0. + */ + +#include +#include +#include + +extern double jn_error(); + +double +jn(n, x) +register int n; +register double x; +{ + double a, b, temp, t; + int i; + + if (_ABS(x) > X_TLOSS) + return (jn_error(n, x, 1)); + if (n == 0) + return (j0(x)); + if (x == 0) + return (x); + if (n < 0) { + n = -n; + x = -x; + } + if (n == 1) + return (j1(x)); + if (n <= x) { + a = j0(x); + b = j1(x); + for (i = 1; i < n; i++) { + temp = b; + b = (i + i)/x * b - a; + a = temp; + } + return (b); + } + temp = x * x; + for (t = 0, i = n + 16; i > n; i--) + t = temp/(i + i - t); + a = t = x/(n + n - t); + b = 1; + for (i = n - 1; i > 0; i--) { + temp = b; + b = (i + i)/x * b - a; + a = temp; + } + return (t * j0(x)/b); +} + +double +yn(n, x) +register int n; +register double x; +{ + double a, b, temp; + int i, neg; + + if (x <= 0) { + struct exception exc; + + exc.type = DOMAIN; + exc.name = "yn"; + exc.arg1 = n; + exc.arg2 = x; + exc.retval = -HUGE; + if (!matherr(&exc)) { + (void) write(2, "yn: DOMAIN error\n", 17); + errno = EDOM; + } + return (exc.retval); + } + if (x > X_TLOSS) + return (jn_error(n, x, 0)); + if (n == 0) + return (y0(x)); + neg = 0; + if (n < 0) { + n = -n; + neg = n % 2; + } + b = y1(x); + if (n > 1) { + a = y0(x); + for (i = 1; i < n; i++) { + temp = b; + b = (i + i)/x * b - a; + a = temp; + } + } + return (neg ? -b : b); +} + +static double +jn_error(n, x, jnflag) +int n; +double x; +int jnflag; +{ + struct exception exc; + + exc.type = TLOSS; + exc.name = jnflag ? "jn" : "yn"; + exc.arg1 = n; + exc.arg2 = x; + exc.retval = 0.0; + if (!matherr(&exc)) { + (void) write(2, exc.name, 2); + (void) write(2, ": TLOSS error\n", 14); + errno = ERANGE; + } + return (exc.retval); +} diff --git a/WEGA/contrib/lcc/lib/libm/log.c b/WEGA/contrib/lcc/lib/libm/log.c new file mode 100644 index 0000000..df6dabc --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/log.c @@ -0,0 +1,88 @@ +/* ~(#)log.c 2.1 */ +/* ~(#)log.c 1.12 */ +/*LINTLIBRARY*/ +/* + * log returns the natural logarithm of its double-precision argument. + * log10 returns the base-10 logarithm of its double-precision argument. + * Returns EDOM error and value -HUGE if argument <= 0. + * Algorithm and coefficients from Cody and Waite (1980). + * Calls frexp. + */ + +#include +#include +#define C1 0.693359375 +#define C2 -2.121944400546905827679e-4 + +extern double log_error(); + +double +log(x) +register double x; +{ + static double p[] = { + -0.78956112887491257267e0, + 0.16383943563021534222e2, + -0.64124943423745581147e2, + }, q[] = { + 1.0, + -0.35667977739034646171e2, + 0.31203222091924532844e3, + -0.76949932108494879777e3, + }; + register double y; + int n; /* can't be in register because of frexp() below */ + + if (x <= 0) + return (log_error(x, "log", 3)); + y = 1.0; + x = frexp(x, &n); + if (x < M_SQRT1_2) { + n--; + y = 0.5; + } + x = (x - y)/(x + y); + x += x; + y = x * x; + x += x * y * _POLY2(y, p)/_POLY3(y, q); + y = (double)n; + x += y * C2; + return (x + y * C1); +} + +double +log10(x) +register double x; +{ + return (x > 0 ? log(x) * M_LOG10E : log_error(x, "log10", 5)); +} + +static double +log_error(x, f_name, name_len) +double x; +char *f_name; +unsigned int name_len; +{ + register char *err_type; + unsigned int mess_len; + struct exception exc; + + exc.name = f_name; + exc.retval = -HUGE; + exc.arg1 = x; + if (x) { + exc.type = DOMAIN; + err_type = ": DOMAIN error\n"; + mess_len = 15; + } else { + exc.type = SING; + err_type = ": SING error\n"; + mess_len = 13; + } + if (!matherr(&exc)) { + (void) write(2, f_name, name_len); + (void) write(2, err_type, mess_len); + errno = EDOM; + } + return (exc.retval); +} diff --git a/WEGA/contrib/lcc/lib/libm/matherr.c b/WEGA/contrib/lcc/lib/libm/matherr.c new file mode 100644 index 0000000..8fd181b --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/matherr.c @@ -0,0 +1,13 @@ +/* ~(#)matherr.c 2.1 */ +/* ~(#)matherr.c 1.2 */ +/*LINTLIBRARY*/ + +#include + +/*ARGSUSED*/ +int +matherr(x) +struct exception *x; +{ + return (0); +} diff --git a/WEGA/contrib/lcc/lib/libm/pow.c b/WEGA/contrib/lcc/lib/libm/pow.c new file mode 100644 index 0000000..e1b96d6 --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/pow.c @@ -0,0 +1,89 @@ +/* ~(#)pow.c 2.1 */ +/* ~(#)pow.c 1.14 */ +/*LINTLIBRARY*/ +/* + * pow(x, y) returns x ** y. + * Returns EDOM error and value 0 for 0 to a non-positive power + * or negative to a non-integral power; + * ERANGE error and value HUGE or -HUGE when the correct value + * would overflow, or 0 when the correct value would underflow. + * uses log and exp + * This version accepts negative x and integral y, + * apparently in violation of the ANSI FORTRAN standard for x ** y. + * There is a much more accurate but much more elaborate algorithm + * in Cody and Waite, which should be substituted for this. + */ + +#include +#include +#include + +double +pow(x, y) +register double x, y; +{ + register int neg; + struct exception exc; + + if (y == 1) /* easy special case */ + return (x); + exc.name = "pow"; + exc.arg1 = x; + exc.arg2 = y; + exc.retval = 0.0; + if (!x) { + if (y > 0) + return (x); /* (0 ** pos) == 0 */ + goto domain; + } + neg = 0; + if (x < 0) { /* test using integer arithmetic if possible */ + if (y >= -MAXLONG && y <= MAXLONG) { + register long ly = (long)y; + + if ((double)ly != y) + goto domain; /* y not integral */ + neg = ly % 2; + } else { + double fr, dum, modf(); + + if (fr = modf(0.5 * y, &dum)) { + if (fr != 0.5) + goto domain; /* y not integral */ + neg++; /* y is odd */ + } + } + x = -x; + } + if (x != 1) { /* x isn't the final result */ + /* the following code protects against multiplying x and y + * until there is no chance of multiplicative overflow */ + if ((x = log(x)) < 0) { /* preserve sign of product */ + x = -x; + y = -y; + } + if (y > LN_MAXDOUBLE/x) { + exc.type = OVERFLOW; + exc.retval = neg ? -HUGE : HUGE; + if (!matherr(&exc)) + errno = ERANGE; + return (exc.retval); + } + if (y < LN_MINDOUBLE/x) { + exc.type = UNDERFLOW; + if (!matherr(&exc)) + errno = ERANGE; + return (exc.retval); + } + x = exp(x * y); /* finally; no mishap can occur */ + } + return (neg ? -x : x); + +domain: + exc.type = DOMAIN; + if (!matherr(&exc)) { + (void) write(2, "pow: DOMAIN error\n", 18); + errno = EDOM; + } + return (exc.retval); +} diff --git a/WEGA/contrib/lcc/lib/libm/sin.c b/WEGA/contrib/lcc/lib/libm/sin.c new file mode 100644 index 0000000..d75a26e --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/sin.c @@ -0,0 +1,106 @@ +/* ~(#)sin.c 2.1 */ +/* ~(#)sin.c 1.13 */ +/*LINTLIBRARY*/ +/* + * sin and cos of double-precision argument. + * Returns ERANGE error and value 0 if argument too large. + * Algorithm and coefficients from Cody and Waite (1980). + */ + +#include +#include +#include + +double +sin(x) +double x; +{ + extern double sin_cos(); + + return (sin_cos(x, 0)); +} + +static double +sin_cos(x, cosflag) +register double x; +int cosflag; +{ + register double y; + register int neg = 0; + struct exception exc; + + exc.arg1 = x; + if (x < 0) { + x = -x; + neg++; + } + y = x; + if (cosflag) { + neg = 0; + y += M_PI_2; + exc.name = "cos"; + } else + exc.name = "sin"; + if (y > X_TLOSS) { + exc.type = TLOSS; + exc.retval = 0.0; + if (!matherr(&exc)) { + (void) write(2, exc.name, 3); + (void) write(2, ": TLOSS error\n", 14); + errno = ERANGE; + } + return (exc.retval); + } + y = y * M_1_PI + 0.5; + if (x <= MAXLONG) { /* reduce using integer arithmetic if possible */ + register long n = (long)y; + + y = (double)n; + if (cosflag) + y -= 0.5; + _REDUCE(long, x, y, 3.1416015625, -8.908910206761537356617e-6); + neg ^= (int)n % 2; + } else { + extern double modf(); + double dn; + + x = (modf(y, &dn) - 0.5) * M_PI; + if (modf(0.5 * dn, &dn)) + neg ^= 1; + } + if (x < 0) { + x = -x; + neg ^= 1; + } + if (x > X_EPS) { /* skip for efficiency and to prevent underflow */ + static double p[] = { + 0.27204790957888846175e-14, + -0.76429178068910467734e-12, + 0.16058936490371589114e-9, + -0.25052106798274584544e-7, + 0.27557319210152756119e-5, + -0.19841269841201840457e-3, + 0.83333333333331650314e-2, + -0.16666666666666665052e0, + }; + + y = x * x; + x += x * y * _POLY7(y, p); + if (x > 1.0) /* inhibit roundoff out of range */ + x = 1.0; + } + exc.retval = neg ? -x : x; + if (exc.arg1 > X_PLOSS || exc.arg1 < -X_PLOSS) { + exc.type = PLOSS; + if (!matherr(&exc)) + errno = ERANGE; + } + return (exc.retval); +} + +double +cos(x) +register double x; +{ + return (x ? sin_cos(x, 1) : 1.0); +} diff --git a/WEGA/contrib/lcc/lib/libm/sinh.c b/WEGA/contrib/lcc/lib/libm/sinh.c new file mode 100644 index 0000000..70ee8f4 --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/sinh.c @@ -0,0 +1,91 @@ +/* ~(#)sinh.c 2.1 */ +/* ~(#)sinh.c 1.9 */ +/*LINTLIBRARY*/ +/* + * sinh returns the hyperbolic sine of its double-precision argument. + * A series is used for arguments smaller in magnitude than 1. + * The exponential function is used for arguments + * greater in magnitude than 1. + * + * cosh returns the hyperbolic cosine of its double-precision argument. + * cosh is computed from the exponential function for + * all arguments. + * + * Returns ERANGE error and value HUGE (or -HUGE for sinh of + * negative argument) if the correct result would overflow. + * + * Algorithm and coefficients from Cody and Waite (1980). + */ + +#include +#include +#include +#define X_MAX (LN_MAXDOUBLE + M_LN2) +#define LNV 0.6931610107421875 +#define V2M1 0.13830277879601902638e-4 + +extern double sinh_exc(); + +double +sinh(x) +register double x; +{ + register double y = _ABS(x); + + if (y <= 1) { + static double p[] = { + -0.78966127417357099479e0, + -0.16375798202630751372e3, + -0.11563521196851768270e5, + -0.35181283430177117881e6, + }, q[] = { + 1.0, + -0.27773523119650701667e3, + 0.36162723109421836460e5, + -0.21108770058106271242e7, + }; + + if (y < X_EPS) /* for efficiency and to prevent underflow */ + return (x); + y = x * x; + return (x + x * y * _POLY3(y, p)/_POLY3(y, q)); + } + if (y > LN_MAXDOUBLE) /* exp(x) would overflow */ + return (sinh_exc(x, y, 1)); + x = exp(x); + return (0.5 * (x - 1.0/x)); +} + +double +cosh(x) +register double x; +{ + register double y = _ABS(x); + + if (y > LN_MAXDOUBLE) /* exp(x) would overflow */ + return (sinh_exc(x, y, 0)); + x = exp(y); + return (0.5 * (x + 1.0/x)); +} + +static double +sinh_exc(x, y, sinhflag) +register double x, y; +register int sinhflag; +{ + int neg = (x < 0 && sinhflag); /* sinh of negative argument */ + struct exception exc; + + if (y < X_MAX) { /* result is still representable */ + x = exp(y - LNV); + x += V2M1 * x; + return (neg ? -x : x); + } + exc.type = OVERFLOW; + exc.name = sinhflag ? "sinh" : "cosh"; + exc.arg1 = x; + exc.retval = neg ? -HUGE : HUGE; + if (!matherr(&exc)) + errno = ERANGE; + return (exc.retval); +} diff --git a/WEGA/contrib/lcc/lib/libm/sqrt.c b/WEGA/contrib/lcc/lib/libm/sqrt.c new file mode 100644 index 0000000..88835c6 --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/sqrt.c @@ -0,0 +1,48 @@ +/* ~(#)sqrt.c 2.1 */ +/* ~(#)sqrt.c 1.9 */ +/*LINTLIBRARY*/ +/* + * sqrt returns the square root of its double-precision argument, + * using Newton's method. + * Returns EDOM error and value 0 if argument negative. + * Calls frexp and ldexp. + */ + +#include +#include +#define ITERATIONS 4 + +double +sqrt(x) +register double x; +{ + register double y; + int iexp; /* can't be in register because of frexp() below */ + register int i = ITERATIONS; + + if (x <= 0) { + struct exception exc; + + if (!x) + return (x); /* sqrt(0) == 0 */ + exc.type = DOMAIN; + exc.name = "sqrt"; + exc.arg1 = x; + exc.retval = 0.0; + if (!matherr(&exc)) { + (void) write(2, "sqrt: DOMAIN error\n", 19); + errno = EDOM; + } + return (exc.retval); + } + y = frexp(x, &iexp); /* 0.5 <= y < 1 */ + if (iexp % 2) { /* iexp is odd */ + --iexp; + y += y; /* 1 <= y < 2 */ + } + y = ldexp(y + 1.0, iexp/2 - 1); /* first guess for sqrt */ + do { + y = 0.5 * (y + x/y); + } while (--i > 0); + return (y); +} diff --git a/WEGA/contrib/lcc/lib/libm/tan.c b/WEGA/contrib/lcc/lib/libm/tan.c new file mode 100644 index 0000000..41e9821 --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/tan.c @@ -0,0 +1,90 @@ +/* ~(#)tan.c 2.1 */ +/* ~(#)tan.c 1.12 */ +/*LINTLIBRARY*/ +/* + * tan returns the tangent of its double-precision argument. + * Returns ERANGE error and value 0 if argument too large. + * Algorithm and coefficients from Cody and Waite (1980). + */ + +#include +#include +#include + +double +tan(x) +register double x; +{ + register double y; + register int neg = 0, invert; + struct exception exc; + + exc.name = "tan"; + exc.arg1 = x; + if (x < 0) { + x = -x; + neg++; + } + if (x > X_TLOSS/2) { + exc.type = TLOSS; + exc.retval = 0.0; + if (!matherr(&exc)) { + (void) write(2, "tan: TLOSS error\n", 17); + errno = ERANGE; + } + return (exc.retval); + } + y = x * M_2_PI + 0.5; + if (x <= MAXLONG) { /* reduce using integer arithmetic if possible */ + register long n = (long)y; + + y = (double)n; + _REDUCE(long, x, y, 1.57080078125, -4.454455103380768678308e-6); + invert = (int)n % 2; + } else { + extern double modf(); + double dn; + + x = (modf(y, &dn) - 0.5) * M_PI_2; + invert = modf(0.5 * dn, &dn) ? 1 : 0; + } + if (x > -X_EPS && x < X_EPS) { + if (!x && invert) + /* + * This case represents a range reduction + * from (2n + 1) * PI/2 which happens to return 0. + * This is mathematically insignficant, and extremely + * unlikely in any case, so it is best simply + * to substitute an arbitrarily small value + * as if the reduction were one bit away from 0, + * and not to deal with this as a range error. + */ + x = 1.0/X_TLOSS; + y = 1.0; /* skip for efficiency and to prevent underflow */ + } else { + static double p[] = { + -0.17861707342254426711e-4, + 0.34248878235890589960e-2, + -0.13338350006421960681e+0, + }, q[] = { + 0.49819433993786512270e-6, + -0.31181531907010027307e-3, + 0.25663832289440112864e-1, + -0.46671683339755294240e+0, + 1.0, + }; + + y = x * x; + x += x * y * _POLY2(y, p); + y = _POLY4(y, q); + } + if (neg) + x = -x; + exc.retval = invert ? -y/x : x/y; + if (exc.arg1 > X_PLOSS/2 || exc.arg1 < -X_PLOSS/2) { + exc.type = PLOSS; + if (!matherr(&exc)) + errno = ERANGE; + } + return (exc.retval); +} diff --git a/WEGA/contrib/lcc/lib/libm/tanh.c b/WEGA/contrib/lcc/lib/libm/tanh.c new file mode 100644 index 0000000..f6181de --- /dev/null +++ b/WEGA/contrib/lcc/lib/libm/tanh.c @@ -0,0 +1,47 @@ +/* ~(#)tanh.c 2.1 */ +/* ~(#)tanh.c 1.7 */ +/*LINTLIBRARY*/ +/* + * tanh returns the hyperbolic tangent of its double-precision argument. + * It calls exp for absolute values of the argument > ~0.55. + * There are no error returns. + * Algorithm and coefficients from Cody and Waite (1980). + */ + +#include +#include +#define X_BIG (0.5 * (LN_MAXDOUBLE + M_LN2)) +#define LN_3_2 0.54930614433405484570 + +double +tanh(x) +register double x; +{ + register int neg = 0; + + if (x < 0) { + x = -x; + neg++; + } + if (x > X_BIG) + x = 1.0; + else if (x > LN_3_2) { + x = 0.5 - 1.0/(exp(x + x) + 1.0); /* two steps recommended */ + x += x; /* for wobbling-precision machines (like IBM) */ + } else if (x > X_EPS) { /* skip for efficiency and to prevent underflow */ + static double p[] = { + -0.96437492777225469787e0, + -0.99225929672236083313e2, + -0.16134119023996228053e4, + }, q[] = { + 1.0, + 0.11274474380534949335e3, + 0.22337720718962312926e4, + 0.48402357071988688686e4, + }; + register double y = x * x; + + x += x * y * _POLY2(y, p)/_POLY3(y, q); + } + return (neg ? -x : x); +}