From: ceriel Date: Mon, 25 Jul 1988 11:26:26 +0000 (+0000) Subject: replaced mathematical routines by our own X-Git-Tag: release-5-5~3003 X-Git-Url: https://git.ndcode.org/public/gitweb.cgi?a=commitdiff_plain;h=324c95ae6260e0edcf2de6833ce38e45f2f6986b;p=ack.git replaced mathematical routines by our own --- diff --git a/lang/pc/libpc/atn.c b/lang/pc/libpc/atn.c index 2f88b9780..1358e4db2 100644 --- a/lang/pc/libpc/atn.c +++ b/lang/pc/libpc/atn.c @@ -1,74 +1,102 @@ -/* $Header$ */ - /* - floating-point arctangent + * (c) copyright 1988 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 + */ - atan returns the value of the arctangent of its - argument in the range [-pi/2,pi/2]. - - there are no error returns. +/* $Header$ */ - coefficients are #5077 from Hart & Cheney. (19.56D) -*/ +#include +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) -static double sq2p1 = 2.414213562373095048802e0; -static double sq2m1 = .414213562373095048802e0; -static double pio2 = 1.570796326794896619231e0; -static double pio4 = .785398163397448309615e0; -static double p4 = .161536412982230228262e2; -static double p3 = .26842548195503973794141e3; -static double p2 = .11530293515404850115428136e4; -static double p1 = .178040631643319697105464587e4; -static double p0 = .89678597403663861959987488e3; -static double q4 = .5895697050844462222791e2; -static double q3 = .536265374031215315104235e3; -static double q2 = .16667838148816337184521798e4; -static double q1 = .207933497444540981287275926e4; -static double q0 = .89678597403663861962481162e3; + arctan(x) = arctan(xi) + arctan(t) -/* - xatan evaluates a series valid in the - range [-0.414...,+0.414...]. -*/ + For the interval [0, p/16] an approximation is used: + arctan(x) = x * P(x*x)/Q(x*x) + */ + 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 }}; -static double -xatan(arg) -double arg; -{ - double argsq; - double value; + /* Hart & Cheney # 5037 */ - argsq = arg*arg; - value = ((((p4*argsq + p3)*argsq + p2)*argsq + p1)*argsq + p0); - value = value/(((((argsq + q4)*argsq + q3)*argsq + q2)*argsq + q1)*argsq + q0); - return(value*arg); -} + static double p[5] = { + 0.7698297257888171026986294745e+03, + 0.1557282793158363491416585283e+04, + 0.1033384651675161628243434662e+04, + 0.2485841954911840502660889866e+03, + 0.1566564964979791769948970100e+02 + }; -static double -satan(arg) -double arg; -{ - if(arg < sq2m1) - return(xatan(arg)); - else if(arg > sq2p1) - return(pio2 - xatan(1/arg)); - else - return(pio4 + xatan((arg-1)/(arg+1))); -} + static double q[6] = { + 0.7698297257888171026986294911e+03, + 0.1813892701754635858982709369e+04, + 0.1484049607102276827437401170e+04, + 0.4904645326203706217748848797e+03, + 0.5593479839280348664778328000e+02, + 0.1000000000000000000000000000e+01 + }; + int negative = x < 0.0; + register struct precomputed *pr = prec; -/* - atan makes its argument positive and - calls the inner routine satan. -*/ + if (negative) { + 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)); + } + else { + double xsq = x*x; -double -_atn(arg) -double arg; -{ - if(arg>0) - return(satan(arg)); - else - return(-satan(-arg)); + x = x * POLYNOM4(xsq, p)/POLYNOM5(xsq, q); + } + return negative ? -x : x; } diff --git a/lang/pc/libpc/exp.c b/lang/pc/libpc/exp.c index 8181f1a14..3ad52e640 100644 --- a/lang/pc/libpc/exp.c +++ b/lang/pc/libpc/exp.c @@ -1,106 +1,112 @@ -/* $Header$ */ - -#include - -extern double _fif(); -extern double _fef(); -extern _trp(); - /* - exp returns the exponential function of its - floating-point argument. - - The coefficients are #1069 from Hart and Cheney. (22.35D) -*/ + * (c) copyright 1988 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 + */ -#define HUGE 1.701411733192644270e38 +/* $Header$ */ -static double p0 = .2080384346694663001443843411e7; -static double p1 = .3028697169744036299076048876e5; -static double p2 = .6061485330061080841615584556e2; -static double q0 = .6002720360238832528230907598e7; -static double q1 = .3277251518082914423057964422e6; -static double q2 = .1749287689093076403844945335e4; -static double log2e = 1.4426950408889634073599247; -static double sqrt2 = 1.4142135623730950488016887; -static double maxf = 10000.0; +#include +#include +extern _trp(); static double -floor(d) -double d; +floor(x) + double x; { - if (d<0) { - d = -d; - if (_fif(d, 1.0, &d) != 0) - d += 1; - d = -d; - } else - _fif(d, 1.0, &d); - return(d); + 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(fr,exp) -double fr; -int exp; +ldexp(fl,exp) + double fl; + int exp; { - int neg,i; + extern double _fef(); + int sign = 1; + int currexp; - neg = 1; - if (fr < 0) { - fr = -fr; - neg = -1; - } - fr = _fef(fr, &i); - /* - while (fr < 0.5) { - fr *= 2; - exp--; + if (fl<0) { + fl = -fl; + sign = -1; } - */ - exp += i; - if (exp > 127) { - _trp(EEXP); - return(neg * HUGE); - } - if (exp < -127) - return(0); - while (exp > 14) { - fr *= (1<<14); - exp -= 14; + fl = _fef(fl,&currexp); + exp += currexp; + if (exp > 0) { + while (exp>30) { + fl *= (double) (1L << 30); + exp -= 30; + } + fl *= (double) (1L << exp); } - while (exp < -14) { - fr /= (1<<14); - exp += 14; + else { + while (exp<-30) { + fl /= (double) (1L << 30); + exp += 30; + } + fl /= (double) (1L << -exp); } - if (exp > 0) - fr *= (1< maxf) { - _trp(EEXP); - return(HUGE); + if (x < M_LN_MIN_D) { + return M_MIN_D; + } + if (x >= M_LN_MAX_D) { + if (x > M_LN_MAX_D) { + _trp(EEXP); + return HUGE; + } + 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; } - arg *= log2e; - ent = floor(arg); - fract = (arg-ent) - 0.5; - xsq = fract*fract; - temp1 = ((p2*xsq+p1)*xsq+p0)*fract; - temp2 = ((xsq+q2)*xsq+q1)*xsq + q0; - return(ldexp(sqrt2*(temp2+temp1)/(temp2-temp1), ent)); + 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; } diff --git a/lang/pc/libpc/log.c b/lang/pc/libpc/log.c index 71aa2e6f1..c3db6e910 100644 --- a/lang/pc/libpc/log.c +++ b/lang/pc/libpc/log.c @@ -1,59 +1,55 @@ -/* $Header$ */ - -#include - -extern double _fef(); -extern _trp(); - /* - log returns the natural logarithm of its floating - point argument. + * (c) copyright 1988 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 + */ - The coefficients are #2705 from Hart & Cheney. (19.38D) - - It calls _fef. -*/ - -#define HUGE 1.701411733192644270e38 +/* $Header$ */ -static double log2 = 0.693147180559945309e0; -static double sqrto2 = 0.707106781186547524e0; -static double p0 = -.240139179559210510e2; -static double p1 = 0.309572928215376501e2; -static double p2 = -.963769093368686593e1; -static double p3 = 0.421087371217979714e0; -static double q0 = -.120069589779605255e2; -static double q1 = 0.194809660700889731e2; -static double q2 = -.891110902798312337e1; +#include +#include +extern _trp(); double -_log(arg) -double arg; +_log(x) + double x; { - double x,z, zsq, temp; - int exp; - - if(arg <= 0) { - _trp(ELOG); - return(-HUGE); - } - x = _fef(arg,&exp); - /* - while(x < 0.5) { - x =* 2; - exp--; - } + /* log(x) = z*P(z*z)/Q(z*z), z = (x-1)/(x+1), x in [1/sqrt(2), sqrt(2)] */ - if(x static double -sinus(arg, quad) -double arg; -int quad; +sinus(x, quadrant) + double x; { - double e, f; - double ysq; - double x,y; - int k; - double temp1, temp2; + /* 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 + }; - x = arg; - if(x<0) { + 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; + + if (x < 0) { + quadrant += 2; x = -x; - quad = quad + 2; } - x = x*twoopi; /*underflow?*/ - if(x>32764){ - y = _fif(x, 10.0, &e); - e = e + quad; - _fif(0.25, e, &f); - quad = e - 4*f; - }else{ - k = x; - y = x - k; - quad = (quad + k) & 03; + 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 (quad & 01) - y = 1-y; - if(quad > 1) - y = -y; + 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; + x += x2; + x -= n * A2; +#undef A1 +#undef A2 + } + else { + extern double _fif(); + double dummy; - ysq = y*y; - temp1 = ((((p4*ysq+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y; - temp2 = ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0); - return(temp1/temp2); + 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) { + x = -x; + } + xsqr = x * x; + x = x * POLYNOM5(xsqr, p) / POLYNOM5(xsqr, q); + return x; } double -_cos(arg) -double arg; +_sin(x) + double x; { - if(arg<0) - arg = -arg; - return(sinus(arg, 1)); + return sinus(x, 0); } double -_sin(arg) -double arg; +_cos(x) + double x; { - return(sinus(arg, 0)); + if (x < 0) x = -x; + return sinus(x, 1); } diff --git a/lang/pc/libpc/sqt.c b/lang/pc/libpc/sqt.c index 49803aa27..39c06cfcd 100644 --- a/lang/pc/libpc/sqt.c +++ b/lang/pc/libpc/sqt.c @@ -1,60 +1,72 @@ +/* + * (c) copyright 1988 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 + */ + /* $Header$ */ -#include +#include +#include +extern _trp(); -extern double _fef(); -extern _trp(); +#define NITER 5 -/* - sqrt returns the square root of its floating - point argument. Newton's method. +static double +ldexp(fl,exp) + double fl; + int exp; +{ + extern double _fef(); + int sign = 1; + int currexp; - calls _fef -*/ + 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 -_sqt(arg) -double arg; +_sqt(x) + double x; { - double x, temp; - int exp; - int i; + extern double _fef(); + int exponent; + double val; - if(arg <= 0) { - if(arg < 0) - _trp(ESQT); - return(0); - } - x = _fef(arg,&exp); - /* - while(x < 0.5) { - x =* 2; - exp--; - } - */ - /* - * NOTE - * this wont work on 1's comp - */ - if(exp & 1) { - x *= 2; - exp--; + if (x <= 0) { + if (x < 0) _trp(ESQT); + return 0; } - temp = 0.5*(1 + x); - while(exp > 28) { - temp *= (1<<14); - exp -= 28; + val = _fef(x, &exponent); + if (exponent & 1) { + exponent--; + val *= 2; } - while(exp < -28) { - temp /= (1<<14); - exp += 28; + val = ldexp(val + 1.0, exponent/2 - 1); + /* was: val = (val + 1.0)/2.0; val = ldexp(val, exponent/2); */ + for (exponent = NITER - 1; exponent >= 0; exponent--) { + val = (val + x / val) / 2.0; } - if(exp >= 0) - temp *= 1 << (exp/2); - else - temp /= 1 << (-exp/2); - for(i=0; i<=4; i++) - temp = 0.5*(temp + arg/temp); - return(temp); + return val; }