diff --git a/lang/pc/libpc/atn.c b/lang/pc/libpc/atn.c index eb09b9283..92205a54d 100644 --- a/lang/pc/libpc/atn.c +++ b/lang/pc/libpc/atn.c @@ -6,6 +6,7 @@ */ /* $Header$ */ + #define __NO_DEFS #include @@ -13,90 +14,55 @@ double _atn(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[] = { + 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 */ }; - static double q[6] = { - 0.7698297257888171026986294911e+03, - 0.1813892701754635858982709369e+04, - 0.1484049607102276827437401170e+04, - 0.4904645326203706217748848797e+03, - 0.5593479839280348664778328000e+02, - 0.1000000000000000000000000000e+01 - }; + int neg = x < 0; + int n; + double g; - int negative = x < 0.0; - register struct precomputed *pr = prec; - - if (negative) { + if (neg) { x = -x; } - while (x > pr->X) pr++; - if (pr != prec) { - x = pr->arctan + - _atn(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/pc/libpc/exp.c b/lang/pc/libpc/exp.c index 504b4a13e..f03aacbe4 100644 --- a/lang/pc/libpc/exp.c +++ b/lang/pc/libpc/exp.c @@ -11,76 +11,32 @@ #include extern _trp(); -static double -floor(x) - double x; -{ - extern double _fif(); - double val; - - return _fif(x, 1.0, &val) < 0 ? val - 1.0 : val ; - /* this also works if _fif always returns a positive - fractional part - */ -} - -static double -ldexp(fl,exp) - double fl; - int exp; -{ - extern double _fef(); - int sign = 1; - int currexp; - - if (fl<0) { - fl = -fl; - sign = -1; - } - fl = _fef(fl,&currexp); - exp += currexp; - if (exp > 0) { - while (exp>30) { - fl *= (double) (1L << 30); - exp -= 30; - } - fl *= (double) (1L << exp); - } - else { - while (exp<-30) { - fl /= (double) (1L << 30); - exp += 30; - } - fl /= (double) (1L << -exp); - } - return sign * fl; -} - double _exp(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 }; + double xn, g; + int n; + int negative = x < 0; - int negative = x < 0; - int ipart, large = 0; - double xsqr, xPxx, Qxx; - - if (x < M_LN_MIN_D) { + if (x <= M_LN_MIN_D) { return M_MIN_D; } if (x >= M_LN_MAX_D) { @@ -90,23 +46,24 @@ _exp(x) } return M_MAX_D; } + if (negative) x = -x; + /* ??? 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); + } if (negative) { - x = -x; + g = -g; + n = -n; } - x /= M_LN2; - ipart = floor(x); - x -= ipart; - if (x > 0.5) { - large = 1; - x -= 0.5; - } - 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/pc/libpc/log.c b/lang/pc/libpc/log.c index bde00a299..324fa2103 100644 --- a/lang/pc/libpc/log.c +++ b/lang/pc/libpc/log.c @@ -6,38 +6,34 @@ */ /* $Header$ */ + #define __NO_DEFS #include #include -extern _trp(); 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 b[] = { + -0.76949932108494879777e3, + 0.31203222091924532844e3, + -0.35667977739034646171e2, + 1.0 }; - static double q[5] = { - 0.3752047495388561108727775374e+02, - -0.7979028073715004879439951583e+02, - 0.5616126132118257292058560360e+02, - -0.1450868091858082685362325000e+02, - 0.1000000000000000000000000000e+01 - }; - - extern double _fef(); - double z, zsqr; - int exponent; + extern double _fef(); + double znum, zden, z, w; + int exponent; if (x <= 0) { _trp(ELOG); @@ -45,11 +41,18 @@ _log(x) } x = _fef(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/pc/libpc/sin.c b/lang/pc/libpc/sin.c index cb8a4e030..521b5893d 100644 --- a/lang/pc/libpc/sin.c +++ b/lang/pc/libpc/sin.c @@ -11,90 +11,74 @@ #include 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 */ + /* Algorithm and coefficients from: + "Software manual for the elementary functions" + by W.J. Cody and W. Waite, Prentice-Hall, 1980 + */ - static double p[6] = { - 0.4857791909822798473837058825e+10, - -0.1808816670894030772075877725e+10, - 0.1724314784722489597789244188e+09, - -0.6351331748520454245913645971e+07, - 0.1002087631419532326179108883e+06, - -0.5830988897678192576148973679e+03 + 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 }; - static double q[6] = { - 0.3092566379840468199410228418e+10, - 0.1202384907680254190870913060e+09, - 0.2321427631602460953669856368e+07, - 0.2848331644063908832127222835e+05, - 0.2287602116741682420054505174e+03, - 0.1000000000000000000000000000e+01 - }; - - 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; + extern double _fif(); + + _fif(y, 1.0, &y); + if (_fif(y, 0.5, &x1)) neg = !neg; + if (cos_flag) y -= 0.5; + x2 = _fif(x, 1.0, &x1); + x = x1 - y * A1; x += x2; - x -= n * A2; + x -= y * A2; #undef A1 #undef A2 - } - else { - extern double _fif(); - double dummy; + } - x = _fif(x/M_2PI, 1.0, &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