From 9f7ee118f7ebc41e8519fbbcf1905d65785485c5 Mon Sep 17 00:00:00 2001 From: ceriel Date: Thu, 18 May 1989 15:37:54 +0000 Subject: [PATCH] new versions, mostly from Cody and Waite --- lang/cem/libcc/math/LIST | 1 - lang/cem/libcc/math/asin.c | 54 ++++++++++---- lang/cem/libcc/math/atan.c | 117 +++++++++++-------------------- lang/cem/libcc/math/exp.c | 69 +++++++++--------- lang/cem/libcc/math/log.c | 57 ++++++++------- lang/cem/libcc/math/pow.c | 29 ++++++-- lang/cem/libcc/math/sin.c | 122 ++++++++++++++------------------ lang/cem/libcc/math/sinh.c | 82 ++++++++++++++++------ lang/cem/libcc/math/tan.c | 139 ++++++++++++------------------------- lang/cem/libcc/math/tanh.c | 45 +++++++++--- 10 files changed, 363 insertions(+), 352 deletions(-) diff --git a/lang/cem/libcc/math/LIST b/lang/cem/libcc/math/LIST index 7e0b35ab8..a08ed2910 100644 --- a/lang/cem/libcc/math/LIST +++ b/lang/cem/libcc/math/LIST @@ -3,7 +3,6 @@ asin.c atan2.c atan.c ceil.c -cosh.c fabs.c gamma.c hypot.c diff --git a/lang/cem/libcc/math/asin.c b/lang/cem/libcc/math/asin.c index 0e4b3f05f..6db5cbd9f 100644 --- a/lang/cem/libcc/math/asin.c +++ b/lang/cem/libcc/math/asin.c @@ -12,30 +12,58 @@ extern int errno; - static double asin_acos(x, cosfl) double x; { - int negative = x < 0; - extern double sqrt(), atan(); + int negative = x < 0; + int i; + double g; + extern double sqrt(); + static double p[] = { + -0.27368494524164255994e+2, + 0.57208227877891731407e+2, + -0.39688862997540877339e+2, + 0.10152522233806463645e+2, + -0.69674573447350646411e+0 + }; + static double q[] = { + -0.16421096714498560795e+3, + 0.41714430248260412556e+3, + -0.38186303361750149284e+3, + 0.15095270841030604719e+3, + -0.23823859153670238830e+2, + 1.0 + }; if (negative) { x = -x; } - if (x > 1) { - errno = EDOM; - return 0; + if (x > 0.5) { + i = 1 - cosfl; + if (x > 1) { + errno = EDOM; + return 0; + } + g = 0.5 - 0.5 * y; + y = - sqrt(g); + y += y; } - if (x == 1) { - x = M_PI_2; + else { + /* ??? avoid underflow ??? */ + g = y * y; } - else x = atan(x/sqrt(1-x*x)); - if (negative) x = -x; - if (cosfl) { - return M_PI_2 - x; + y += y * g * POLYNOM4(g, x) / POLYNOM5(g, y); + if (i == 1) { + if (cosfl == 0 || ! negative) { + y = (y + M_PI_4) + M_PI_4; + } + else if (cosfl && negative) { + y = (y + M_PI_2) + M_PI_2; + } } - return x; + if (! cosfl && negative) y = -y; + return y; } double diff --git a/lang/cem/libcc/math/atan.c b/lang/cem/libcc/math/atan.c index 787d9ffbe..15faa580b 100644 --- a/lang/cem/libcc/math/atan.c +++ b/lang/cem/libcc/math/atan.c @@ -10,94 +10,61 @@ #include #include +extern int errno; + double atan(x) double x; { - /* The interval [0, infinity) is treated as follows: - Define partition points Xi - X0 = 0 - X1 = tan(pi/16) - X2 = tan(3pi/16) - X3 = tan(5pi/16) - X4 = tan(7pi/16) - X5 = infinity - and evaluation nodes xi - x2 = tan(2pi/16) - x3 = tan(4pi/16) - x4 = tan(6pi/16) - x5 = infinity - An argument x in [Xn-1, Xn] is now reduced to an argument - t in [-X1, X1] by the following formulas: - - t = 1/xn - (1/(xn*xn) + 1)/((1/xn) + x) - - arctan(x) = arctan(xi) + arctan(t) - - For the interval [0, p/16] an approximation is used: - arctan(x) = x * P(x*x)/Q(x*x) + /* Algorithm and coefficients from: + "Software manual for the elementary functions" + by W.J. Cody and W. Waite, Prentice-Hall, 1980 */ - static struct precomputed { - double X; /* partition point */ - double arctan; /* arctan of evaluation node */ - double one_o_x; /* 1 / xn */ - double one_o_xsq_p_1; /* 1 / (xn*xn) + 1 */ - } prec[5] = { - { 0.19891236737965800691159762264467622, - 0.0, - 0.0, /* these don't matter */ - 0.0 } , - { 0.66817863791929891999775768652308076, /* tan(3pi/16) */ - M_PI_8, - 2.41421356237309504880168872420969808, - 6.82842712474619009760337744841939616 }, - { 1.49660576266548901760113513494247691, /* tan(5pi/16) */ - M_PI_4, - 1.0, - 2.0 }, - { 5.02733949212584810451497507106407238, /* tan(7pi/16) */ - M_3PI_8, - 0.41421356237309504880168872420969808, - 1.17157287525380998659662255158060384 }, - { MAXDOUBLE, - M_PI_2, - 0.0, - 1.0 }}; - - /* Hart & Cheney # 5037 */ - static double p[5] = { - 0.7698297257888171026986294745e+03, - 0.1557282793158363491416585283e+04, - 0.1033384651675161628243434662e+04, - 0.2485841954911840502660889866e+03, - 0.1566564964979791769948970100e+02 + static double p[] = { + -0.13688768894191926929e+2, + -0.20505855195861651981e+2, + -0.84946240351320683534e+1, + -0.83758299368150059274e+0 }; - - static double q[6] = { - 0.7698297257888171026986294911e+03, - 0.1813892701754635858982709369e+04, - 0.1484049607102276827437401170e+04, - 0.4904645326203706217748848797e+03, - 0.5593479839280348664778328000e+02, - 0.1000000000000000000000000000e+01 + static double q[] = { + 0.41066306682575781263e+2, + 0.86157349597130242515e+2, + 0.59578436142597344465e+2, + 0.15024001160028576121e+2, + 1.0 + }; + static double a[] = { + 0.0, + 0.52359877559829887307710723554658381, /* pi/6 */ + M_PI_2, + 1.04719755119659774615421446109316763 /* pi/3 */ }; - int negative = x < 0.0; - register struct precomputed *pr = prec; + int neg = x < 0; + int n; + double g; - if (negative) { + if (neg) { x = -x; } - while (x > pr->X) pr++; - if (pr != prec) { - x = pr->arctan + - atan(pr->one_o_x - pr->one_o_xsq_p_1/(pr->one_o_x + x)); + if (x > 1.0) { + x = 1.0/x; + n = 2; } - else { - double xsq = x*x; + else n = 0; - x = x * POLYNOM4(xsq, p)/POLYNOM5(xsq, q); + if (x > 0.26794919243112270647) { /* 2-sqtr(3) */ + n = n + 1; + x = (((0.73205080756887729353*x-0.5)-0.5)+x)/ + (1.73205080756887729353+x); } - return negative ? -x : x; + + /* ??? avoid underflow ??? */ + + g = x * x; + x += x * g * POLYNOM3(g, p) / POLYNOM4(g, q); + if (n > 1) x = -x; + x += a[n]; + return neg ? -x : x; } diff --git a/lang/cem/libcc/math/exp.c b/lang/cem/libcc/math/exp.c index 2d52ed5c4..b34347ab8 100644 --- a/lang/cem/libcc/math/exp.c +++ b/lang/cem/libcc/math/exp.c @@ -1,5 +1,5 @@ /* - * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + * (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands. * See the copyright notice in the ACK home directory, in the file "Copyright". * * Author: Ceriel J.H. Jacobs @@ -10,32 +10,33 @@ #include #include -extern int errno; +extern int errno; +extern double ldexp(); double exp(x) - double x; + double x; { - /* 2**x = (Q(x*x)+x*P(x*x))/(Q(x*x)-x*P(x*x)) for x in [0,0.5] */ - /* Hart & Cheney #1069 */ + /* Algorithm and coefficients from: + "Software manual for the elementary functions" + by W.J. Cody and W. Waite, Prentice-Hall, 1980 + */ - static double p[3] = { - 0.2080384346694663001443843411e+07, - 0.3028697169744036299076048876e+05, - 0.6061485330061080841615584556e+02 + static double p[] = { + 0.25000000000000000000e+0, + 0.75753180159422776666e-2, + 0.31555192765684646356e-4 }; - static double q[4] = { - 0.6002720360238832528230907598e+07, - 0.3277251518082914423057964422e+06, - 0.1749287689093076403844945335e+04, - 0.1000000000000000000000000000e+01 + static double q[] = { + 0.50000000000000000000e+0, + 0.56817302698551221787e-1, + 0.63121894374398503557e-3, + 0.75104028399870046114e-6 }; - - int negative = x < 0; - int ipart, large = 0; - double xsqr, xPxx, Qxx; - extern double floor(), ldexp(); + double xn, g; + int n; + int negative = x < 0; if (x <= M_LN_MIN_D) { if (x < M_LN_MIN_D) errno = ERANGE; @@ -46,22 +47,18 @@ exp(x) return M_MAX_D; } - if (negative) { - x = -x; - } - x /= M_LN2; - ipart = floor(x); - x -= ipart; - if (x > 0.5) { - large = 1; - x -= 0.5; + /* ??? avoid underflow ??? */ + + n = x * M_LOG2E + 0.5; /* 1/ln(2) = log2(e), 0.5 added for rounding */ + xn = n; + { + double x1 = (long) x; + double x2 = x - x1; + + g = ((x1-xn*0.693359375)+x2) - xn*(-2.1219444005469058277e-4); } - xsqr = x * x; - xPxx = x * POLYNOM2(xsqr, p); - Qxx = POLYNOM3(xsqr, q); - x = (Qxx + xPxx) / (Qxx - xPxx); - if (large) x *= M_SQRT2; - x = ldexp(x, ipart); - if (negative) return 1.0/x; - return x; + xn = g * g; + x = g * POLYNOM2(xn, p); + n += 1; + return (ldexp(0.5 + x/(POLYNOM3(xn, q) - x), n)); } diff --git a/lang/cem/libcc/math/log.c b/lang/cem/libcc/math/log.c index 1b2dc37a2..ac118f230 100644 --- a/lang/cem/libcc/math/log.c +++ b/lang/cem/libcc/math/log.c @@ -1,5 +1,5 @@ /* - * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + * (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands. * See the copyright notice in the ACK home directory, in the file "Copyright". * * Author: Ceriel J.H. Jacobs @@ -10,35 +10,31 @@ #include #include -extern int errno; +extern int errno; +extern double frexp(); double log(x) - double x; + double x; { - /* log(x) = z*P(z*z)/Q(z*z), z = (x-1)/(x+1), x in [1/sqrt(2), sqrt(2)] + /* Algorithm and coefficients from: + "Software manual for the elementary functions" + by W.J. Cody and W. Waite, Prentice-Hall, 1980 */ - /* Hart & Cheney #2707 */ - - static double p[5] = { - 0.7504094990777122217455611007e+02, - -0.1345669115050430235318253537e+03, - 0.7413719213248602512779336470e+02, - -0.1277249755012330819984385000e+02, - 0.3327108381087686938144000000e+00 + static double a[] = { + -0.64124943423745581147e2, + 0.16383943563021534222e2, + -0.78956112887491257267e0 }; - - static double q[5] = { - 0.3752047495388561108727775374e+02, - -0.7979028073715004879439951583e+02, - 0.5616126132118257292058560360e+02, - -0.1450868091858082685362325000e+02, - 0.1000000000000000000000000000e+01 + static double b[] = { + -0.76949932108494879777e3, + 0.31203222091924532844e3, + -0.35667977739034646171e2, + 1.0 }; - extern double frexp(); - double z, zsqr; - int exponent; + double znum, zden, z, w; + int exponent; if (x <= 0) { errno = EDOM; @@ -46,11 +42,18 @@ log(x) } x = frexp(x, &exponent); - while (x < M_1_SQRT2) { - x += x; + if (x > M_1_SQRT2) { + znum = (x - 0.5) - 0.5; + zden = x * 0.5 + 0.5; + } + else { + znum = x - 0.5; + zden = znum * 0.5 + 0.5; exponent--; } - z = (x-1)/(x+1); - zsqr = z*z; - return z * POLYNOM4(zsqr, p) / POLYNOM4(zsqr, q) + exponent * M_LN2; + z = znum/zden; w = z * z; + x = z + z * w * (POLYNOM2(w,a)/POLYNOM3(w,b)); + z = exponent; + x += z * (-2.121944400546905827679e-4); + return x + z * 0.693359375; } diff --git a/lang/cem/libcc/math/pow.c b/lang/cem/libcc/math/pow.c index 893f82c4b..4584b43c4 100644 --- a/lang/cem/libcc/math/pow.c +++ b/lang/cem/libcc/math/pow.c @@ -11,13 +11,19 @@ #include extern int errno; +extern double modf(), exp(), log(); double pow(x,y) double x,y; { + /* Simple version for now. The Cody and Waite book has + a very complicated, much more precise version, but + this version has machine-dependant arrays A1 and A2, + and I don't know yet how to solve this ??? + */ double dummy; - extern double modf(), exp(), log(); + int result_neg = 0; if ((x == 0 && y == 0) || (x < 0 && modf(y, &dummy) != 0)) { @@ -28,13 +34,26 @@ pow(x,y) if (x == 0) return x; if (x < 0) { - double val = exp(log(-x) * y); if (modf(y/2.0, &dummy) != 0) { /* y was odd */ - val = - val; + result_neg = 1; } - return val; + x = -x; + } + x = log(x); + if (x < 0) { + x = -x; + y = -y; + } + if (y > M_LN_MAX_D/x) { + errno = ERANGE; + return 0; + } + if (y < M_LN_MIN_D/x) { + errno = ERANGE; + return 0; } - return exp(log(x) * y); + x = exp(x * y); + return result_neg ? -x : x; } diff --git a/lang/cem/libcc/math/sin.c b/lang/cem/libcc/math/sin.c index 97fbdafb4..c7b335bdd 100644 --- a/lang/cem/libcc/math/sin.c +++ b/lang/cem/libcc/math/sin.c @@ -10,93 +10,77 @@ #include #include -extern int errno; +extern int errno; +extern double modf(); static double -sinus(x, quadrant) +sinus(x, cos_flag) double x; { - /* sin(0.5*pi*x) = x * P(x*x)/Q(x*x) for x in [0,1] */ - /* Hart & Cheney # 3374 */ - - static double p[6] = { - 0.4857791909822798473837058825e+10, - -0.1808816670894030772075877725e+10, - 0.1724314784722489597789244188e+09, - -0.6351331748520454245913645971e+07, - 0.1002087631419532326179108883e+06, - -0.5830988897678192576148973679e+03 - }; + /* Algorithm and coefficients from: + "Software manual for the elementary functions" + by W.J. Cody and W. Waite, Prentice-Hall, 1980 + */ - static double q[6] = { - 0.3092566379840468199410228418e+10, - 0.1202384907680254190870913060e+09, - 0.2321427631602460953669856368e+07, - 0.2848331644063908832127222835e+05, - 0.2287602116741682420054505174e+03, - 0.1000000000000000000000000000e+01 + static double r[] = { + -0.16666666666666665052e+0, + 0.83333333333331650314e-2, + -0.19841269841201840457e-3, + 0.27557319210152756119e-5, + -0.25052106798274584544e-7, + 0.16058936490371589114e-9, + -0.76429178068910467734e-12, + 0.27204790957888846175e-14 }; - double xsqr; - int t; + double xsqr; + double y; + int neg = 0; if (x < 0) { - quadrant += 2; x = -x; + neg = 1; } - if (M_PI_2 - x == M_PI_2) { - switch(quadrant) { - case 0: - case 2: - return 0.0; - case 1: - return 1.0; - case 3: - return -1.0; - } + if (cos_flag) { + neg = 0; + y = M_PI_2 + x; } - if (x >= M_2PI) { - if (x <= 0x7fffffff) { - /* Use extended precision to calculate reduced argument. - Split 2pi in 2 parts a1 and a2, of which the first only - uses some bits of the mantissa, so that n * a1 is - exactly representable, where n is the integer part of - x/pi. - Here we used 12 bits of the mantissa for a1. - Also split x in integer part x1 and fraction part x2. - We then compute x-n*2pi as ((x1 - n*a1) + x2) - n*a2. - */ -#define A1 6.2822265625 -#define A2 0.00095874467958647692528676655900576 - double n = (long) (x / M_2PI); - double x1 = (long) x; - double x2 = x - x1; - x = x1 - n * A1; + else y = x; + + /* ??? avoid loss of significance, if y is too large, error ??? */ + + y = y * M_1_PI + 0.5; + + /* Use extended precision to calculate reduced argument. + Here we used 12 bits of the mantissa for a1. + Also split x in integer part x1 and fraction part x2. + */ +#define A1 3.1416015625 +#define A2 -8.908910206761537356617e-6 + { + double x1, x2; + + modf(y, &y); + if (modf(0.5*y, &x1)) neg = !neg; + if (cos_flag) y -= 0.5; + x2 = modf(x, &x1); + x = x1 - y * A1; x += x2; - x -= n * A2; + x -= y * A2; #undef A1 #undef A2 - } - else { - extern double modf(); - double dummy; - - x = modf(x/M_2PI, &dummy) * M_2PI; - } } - x /= M_PI_2; - t = x; - x -= t; - quadrant = (quadrant + (int)(t % 4)) % 4; - if (quadrant & 01) { - x = 1 - x; - } - if (quadrant > 1) { + + if (x < 0) { + neg = !neg; x = -x; } - xsqr = x * x; - x = x * POLYNOM5(xsqr, p) / POLYNOM5(xsqr, q); - return x; + + /* ??? avoid underflow ??? */ + + y = x * x; + x += x * y * POLYNOM7(y, r); + return neg ? -x : x; } double diff --git a/lang/cem/libcc/math/sinh.c b/lang/cem/libcc/math/sinh.c index f5f94b745..b7e5e7815 100644 --- a/lang/cem/libcc/math/sinh.c +++ b/lang/cem/libcc/math/sinh.c @@ -1,5 +1,5 @@ /* - * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + * (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands. * See the copyright notice in the ACK home directory, in the file "Copyright". * * Author: Ceriel J.H. Jacobs @@ -10,33 +10,73 @@ #include #include -extern int errno; +extern int errno; +extern double exp(); -double -sinh(x) - double x; +static double +sinh_cosh(x, cosh_flag) + double x; { - int negx = x < 0; - extern double exp(); + /* Algorithm and coefficients from: + "Software manual for the elementary functions" + by W.J. Cody and W. Waite, Prentice-Hall, 1980 + */ - if (negx) { - x = -x; + static double p[] = { + -0.35181283430177117881e+6, + -0.11563521196851768270e+5, + -0.16375798202630751372e+3, + -0.78966127417357099479e+0 + }; + static double q[] = { + -0.21108770058106271242e+7, + 0.36162723109421836460e+5, + -0.27773523119650701167e+3, + 1.0 + }; + int negative = x < 0; + double y = negative ? -x : x; + + if (! cosh_flag && y <= 1.0) { + /* ??? check for underflow ??? */ + y = y * y; + return x + x * y * POLYNOM3(y, p)/POLYNOM3(y,q); } - if (x > M_LN_MAX_D) { - /* exp(x) would overflow */ - if (x >= M_LN_MAX_D + M_LN2) { - /* not representable */ - x = HUGE; + + if (y >= M_LN_MAX_D) { + /* exp(y) would cause overflow */ +#define LNV 0.69316101074218750000e+0 +#define VD2M1 0.52820835025874852469e-4 + double w = y - LNV; + + if (w < M_LN_MAX_D+M_LN2-LNV) { + x = exp(w); + x += VD2M1 * x; + } + else { errno = ERANGE; + x = HUGE; } - else x = exp (x - M_LN2); } else { - double expx = exp(x); - x = 0.5 * (expx - 1.0/expx); + double z = exp(y); + + x = 0.5 * (z + (cosh_flag ? 1.0 : -1.0)/z); } - if (negx) { - return -x; - } - return x; + return negative ? -x : x; +} + +double +sinh(x) + double x; +{ + return sinh_cosh(x, 0); +} + +double +cosh(x) + double x; +{ + if (x < 0) x = -x; + return sinh_cosh(x, 1); } diff --git a/lang/cem/libcc/math/tan.c b/lang/cem/libcc/math/tan.c index 350cea46d..36ea0dc1d 100644 --- a/lang/cem/libcc/math/tan.c +++ b/lang/cem/libcc/math/tan.c @@ -10,117 +10,64 @@ #include #include -extern int errno; +extern int errno; +extern double modf(); double tan(x) double x; { - /* First reduce range to [0, pi/4]. - Then use approximation tan(x*pi/4) = x * P(x*x)/Q(x*x). - Hart & Cheney # 4288 - Use: tan(x) = 1/tan(pi/2 - x) - tan(-x) = -tan(x) - tan(x+k*pi) = tan(x) + /* Algorithm and coefficients from: + "Software manual for the elementary functions" + by W.J. Cody and W. Waite, Prentice-Hall, 1980 */ - static double p[5] = { - -0.5712939549476836914932149599e+10, - 0.4946855977542506692946040594e+09, - -0.9429037070546336747758930844e+07, - 0.5282725819868891894772108334e+05, - -0.6983913274721550913090621370e+02 - }; - - static double q[6] = { - -0.7273940551075393257142652672e+10, - 0.2125497341858248436051062591e+10, - -0.8000791217568674135274814656e+08, - 0.8232855955751828560307269007e+06, - -0.2396576810261093558391373322e+04, - 0.1000000000000000000000000000e+01 - }; - int negative = x < 0; - double tmp, tmp1, tmp2; - double xsq; int invert = 0; - int ip; + double y; + static double p[] = { + 1.0, + -0.13338350006421960681e+0, + 0.34248878235890589960e-2, + -0.17861707342254426711e-4 + }; + static double q[] = { + 1.0, + -0.46671683339755294240e+0, + 0.25663832289440112864e-1, + -0.31181531907010027307e-3, + 0.49819433993786512270e-6 + }; if (negative) x = -x; - /* first reduce to [0, pi) */ - if (x >= M_PI) { - if (x <= 0x7fffffff) { - /* Use extended precision to calculate reduced argument. - Split pi in 2 parts a1 and a2, of which the first only - uses some bits of the mantissa, so that n * a1 is - exactly representable, where n is the integer part of - x/pi. - Here we used 12 bits of the mantissa for a1. - Also split x in integer part x1 and fraction part x2. - We then compute x-n*pi as ((x1 - n*a1) + x2) - n*a2. - */ -#define A1 3.14111328125 -#define A2 0.00047937233979323846264338327950288 - double n = (long) (x / M_PI); - double x1 = (long) x; - double x2 = x - x1; - x = x1 - n * A1; + /* ??? avoid loss of significance, error if x is too large ??? */ + + y = x * M_2_PI + 0.5; + + /* Use extended precision to calculate reduced argument. + Here we used 12 bits of the mantissa for a1. + Also split x in integer part x1 and fraction part x2. + */ +#define A1 1.57080078125 +#define A2 -4.454455103380768678308e-6 + { + double x1, x2; + + modf(y, &y); + if (modf(0.5*y, &x1)) invert = 1; + x2 = modf(x, &x1); + x = x1 - y * A1; x += x2; - x -= n * A2; + x -= y * A2; #undef A1 #undef A2 - } - else { - extern double modf(); - - x = modf(x/M_PI, &tmp) * M_PI; - } - } - /* because the approximation uses x*pi/4, we reverse this */ - x /= M_PI_4; - ip = (int) x; - x -= ip; - - switch(ip) { - case 0: - /* [0,pi/4] */ - break; - case 1: - /* [pi/4, pi/2] - tan(x+pi/4) = 1/tan(pi/2 - (x+pi/4)) = 1/tan(pi/4 - x) - */ - invert = 1; - x = 1.0 - x; - break; - case 2: - /* [pi/2, 3pi/4] - tan(x+pi/2) = tan((x+pi/2)-pi) = -tan(pi/2 - x) = - -1/tan(x) - */ - negative = ! negative; - invert = 1; - break; - case 3: - /* [3pi/4, pi) - tan(x+3pi/4) = tan(x-pi/4) = - tan(pi/4-x) - */ - x = 1.0 - x; - negative = ! negative; - break; - } - xsq = x * x; - tmp1 = x*POLYNOM4(xsq, p); - tmp2 = POLYNOM5(xsq, q); - tmp = tmp1 / tmp2; - if (invert) { - if (tmp == 0.0) { - errno = ERANGE; - tmp = HUGE; - } - else tmp = tmp2 / tmp1; } - return negative ? -tmp : tmp; + /* ??? avoid underflow ??? */ + y = x * x; + x += x * y * POLYNOM2(y, p+1); + y = POLYNOM4(y, q); + if (neg) x = -x; + return invert ? -y/x : x/y; } diff --git a/lang/cem/libcc/math/tanh.c b/lang/cem/libcc/math/tanh.c index 3bceb37e9..4da18fb22 100644 --- a/lang/cem/libcc/math/tanh.c +++ b/lang/cem/libcc/math/tanh.c @@ -1,5 +1,5 @@ /* - * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + * (c) copyright 1989 by the Vrije Universiteit, Amsterdam, The Netherlands. * See the copyright notice in the ACK home directory, in the file "Copyright". * * Author: Ceriel J.H. Jacobs @@ -10,18 +10,45 @@ #include #include +extern int errno; +extern double exp(); + double tanh(x) - double x; + double x; { - extern double exp(); + /* Algorithm and coefficients from: + "Software manual for the elementary functions" + by W.J. Cody and W. Waite, Prentice-Hall, 1980 + */ + + static double p[] = { + -0.16134119023996228053e+4, + -0.99225929672236083313e+2, + -0.96437492777225469787e+0 + }; + static double q[] = { + 0.48402357071988688686e+4, + 0.22337720718962312926e+4, + 0.11274474380534949335e+3, + 1.0 + }; + int negative = x < 0; + + if (negative) x = -x; - if (x <= 0.5*M_LN_MIN_D) { - return -1; - } if (x >= 0.5*M_LN_MAX_D) { - return 1; + x = 1.0; + } +#define LN3D2 0.54930614433405484570e+0 /* ln(3)/2 */ + else if (x > LN3D2) { + x = 0.5 - 1.0/(exp(x+x)+1.0); + x += x; + } + else { + /* ??? avoid underflow ??? */ + double g = x*x; + x += x * g * POLYNOM2(g, p)/POLYNOM3(g, q); } - x = exp(x + x); - return (x - 1.0)/(x + 1.0); + return negative ? -x : x; } -- 2.34.1