-/* $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 <math.h>
+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;
}
-/* $Header$ */
-
-#include <pc_err.h>
-
-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 <math.h>
+#include <pc_err.h>
+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<<exp);
- if (exp < 0)
- fr /= (1<<(-exp));
- return(neg * fr);
+ return sign * fl;
}
double
-_exp(arg)
-double arg;
+_exp(x)
+ double x;
{
- double fract;
- double temp1, temp2, xsq;
- int ent;
+ /* 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 */
+
+ static double p[3] = {
+ 0.2080384346694663001443843411e+07,
+ 0.3028697169744036299076048876e+05,
+ 0.6061485330061080841615584556e+02
+ };
+
+ static double q[4] = {
+ 0.6002720360238832528230907598e+07,
+ 0.3277251518082914423057964422e+06,
+ 0.1749287689093076403844945335e+04,
+ 0.1000000000000000000000000000e+01
+ };
+
+ int negative = x < 0;
+ int ipart, large = 0;
+ double xsqr, xPxx, Qxx;
- if(arg == 0)
- return(1);
- if(arg < -maxf)
- return(0);
- if(arg > 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;
}
-/* $Header$ */
-
-#include <pc_err.h>
-
-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 <math.h>
+#include <pc_err.h>
+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<sqrto2) {
- x *= 2;
- exp--;
+ /* Hart & Cheney #2707 */
+
+ static double p[5] = {
+ 0.7504094990777122217455611007e+02,
+ -0.1345669115050430235318253537e+03,
+ 0.7413719213248602512779336470e+02,
+ -0.1277249755012330819984385000e+02,
+ 0.3327108381087686938144000000e+00
+ };
+
+ 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;
+
+ if (x <= 0) {
+ _trp(ELOG);
+ return -HUGE;
}
+ x = _fef(x, &exponent);
+ while (x < M_1_SQRT2) {
+ x += x;
+ exponent--;
+ }
z = (x-1)/(x+1);
- zsq = z*z;
-
- temp = ((p3*zsq + p2)*zsq + p1)*zsq + p0;
- temp = temp/(((zsq + q2)*zsq + q1)*zsq + q0);
- temp = temp*z + exp*log2;
- return(temp);
+ zsqr = z*z;
+ return z * POLYNOM4(zsqr, p) / POLYNOM4(zsqr, q) + exponent * M_LN2;
}
-/* $Header$ */
-
-extern double _fif();
-
/*
- C program for floating point sin/cos.
- Calls _fif.
- There are no error exits.
- Coefficients are #3370 from Hart & Cheney (18.80D).
-*/
+ * (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
+ */
-static double twoopi = 0.63661977236758134308;
-static double p0 = .1357884097877375669092680e8;
-static double p1 = -.4942908100902844161158627e7;
-static double p2 = .4401030535375266501944918e6;
-static double p3 = -.1384727249982452873054457e5;
-static double p4 = .1459688406665768722226959e3;
-static double q0 = .8644558652922534429915149e7;
-static double q1 = .4081792252343299749395779e6;
-static double q2 = .9463096101538208180571257e4;
-static double q3 = .1326534908786136358911494e3;
+/* $Header$ */
+
+#include <math.h>
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);
}
+/*
+ * (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 <pc_err.h>
+#include <math.h>
+#include <pc_err.h>
+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;
}
public / ack.git / commitdiff