--- /dev/null
+/*
+ * (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 <errno.h>
+#include <math.h>
+#include "localmath.h"
+
+
+static double
+asin_acos(double x, int cosfl)
+{
+ int negative = x < 0;
+
+ if (negative) {
+ x = -x;
+ }
+ if (x > 1) {
+ errno = EDOM;
+ return 0;
+ }
+ if (x == 1) {
+ x = M_PI_2;
+ }
+ else x = atan(x/sqrt(1-x*x));
+ if (negative) x = -x;
+ if (cosfl) {
+ return M_PI_2 - x;
+ }
+ return x;
+}
+
+double
+asin(double x)
+{
+ return asin_acos(x, 0);
+}
+
+double
+acos(double x)
+{
+ return asin_acos(x, 1);
+}
--- /dev/null
+/*
+ * (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 <errno.h>
+#include <float.h>
+#include <math.h>
+#include "localmath.h"
+
+double
+atan(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)
+ */
+ 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 },
+ { DBL_MAX,
+ 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 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;
+
+ 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;
+
+ x = x * POLYNOM4(xsq, p)/POLYNOM5(xsq, q);
+ }
+ return negative ? -x : x;
+}
+
--- /dev/null
+/*
+ * (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 <errno.h>
+#include <math.h>
+#include "localmath.h"
+
+double
+atan2(double y, double x)
+{
+ double absx, absy, val;
+
+ if (x == 0 && y == 0) {
+ errno = EDOM;
+ return 0;
+ }
+ absy = y < 0 ? -y : y;
+ absx = x < 0 ? -x : x;
+ if (absy - absx == absy) {
+ /* x negligible compared to y */
+ return y < 0 ? -M_PI_2 : M_PI_2;
+ }
+ if (absx - absy == absx) {
+ /* y negligible compared to x */
+ val = 0.0;
+ }
+ else val = atan(y/x);
+ if (x > 0) {
+ /* first or fourth quadrant; already correct */
+ return val;
+ }
+ if (y < 0) {
+ /* third quadrant */
+ return val - M_PI;
+ }
+ return val + M_PI;
+}
--- /dev/null
+/*
+ * (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 <math.h>
+
+double
+ceil(double x)
+{
+ double val;
+
+ return modf(x, &val) > 0 ? val + 1.0 : val ;
+ /* this also works if modf always returns a positive
+ fractional part
+ */
+}
--- /dev/null
+/*
+ * (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 <errno.h>
+#include <math.h>
+#include "localmath.h"
+
+double
+cosh(double x)
+{
+ if (x < 0) {
+ x = -x;
+ }
+ if (x > M_LN_MAX_D) {
+ /* exp(x) would overflow */
+ if (x >= M_LN_MAX_D + M_LN2) {
+ /* not representable */
+ x = HUGE_VAL;
+ errno = ERANGE;
+ }
+ else x = exp (x - M_LN2);
+ }
+ else {
+ double expx = exp(x);
+ x = 0.5 * (expx + 1.0/expx);
+ }
+ return x;
+}
--- /dev/null
+/*
+ * (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 <errno.h>
+#include <float.h>
+#include <math.h>
+#include "localmath.h"
+
+
+double
+exp(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 */
+
+ 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 (x <= M_LN_MIN_D) {
+ if (x < M_LN_MIN_D) errno = ERANGE;
+ return DBL_MIN;
+ }
+ if (x >= M_LN_MAX_D) {
+ if (x > M_LN_MAX_D) errno = ERANGE;
+ return DBL_MAX;
+ }
+
+ if (negative) {
+ x = -x;
+ }
+ 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;
+}
--- /dev/null
+/*
+ * (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$ */
+
+double
+fabs(double x)
+{
+ return x < 0 ? -x : x;
+}
--- /dev/null
+/*
+ * (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 <math.h>
+
+double
+floor(double x)
+{
+ double val;
+
+ return modf(x, &val) < 0 ? val - 1.0 : val ;
+ /* this also works if modf always returns a positive
+ fractional part
+ */
+}
--- /dev/null
+#
+/*
+ * (c) copyright 1987 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ */
+/* $Header$ */
+
+ mes 2,EM_WSIZE,EM_PSIZE
+#ifndef NOFLOAT
+ exp $frexp
+ pro $frexp,0
+ lal 0
+ loi EM_DSIZE
+ fef EM_DSIZE
+ lal EM_DSIZE
+ loi EM_PSIZE
+ sti EM_WSIZE
+ ret EM_DSIZE
+ end
+#endif
--- /dev/null
+/*
+ * (c) copyright 1987 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ */
+/* $Header$ */
+
+#include <math.h>
+
+double
+ldexp(double fl, int exp)
+{
+ int sign = 1;
+ int currexp;
+
+ if (fl<0) {
+ fl = -fl;
+ sign = -1;
+ }
+ fl = frexp(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;
+}
--- /dev/null
+/*
+ * localmath.h - This header is used by the mathematical library.
+ */
+/* $Header$ */
+
+/* some constants (Hart & Cheney) */
+#define M_PI 3.14159265358979323846264338327950288
+#define M_2PI 6.28318530717958647692528676655900576
+#define M_3PI_4 2.35619449019234492884698253745962716
+#define M_PI_2 1.57079632679489661923132169163975144
+#define M_3PI_8 1.17809724509617246442349126872981358
+#define M_PI_4 0.78539816339744830961566084581987572
+#define M_PI_8 0.39269908169872415480783042290993786
+#define M_1_PI 0.31830988618379067153776752674502872
+#define M_2_PI 0.63661977236758134307553505349005744
+#define M_4_PI 1.27323954473516268615107010698011488
+#define M_E 2.71828182845904523536028747135266250
+#define M_LOG2E 1.44269504088896340735992468100189213
+#define M_LOG10E 0.43429448190325182765112891891660508
+#define M_LN2 0.69314718055994530941723212145817657
+#define M_LN10 2.30258509299404568401799145468436421
+#define M_SQRT2 1.41421356237309504880168872420969808
+#define M_1_SQRT2 0.70710678118654752440084436210484904
+#define M_EULER 0.57721566490153286060651209008240243
+
+/* macros for constructing polynomials */
+#define POLYNOM1(x, a) ((a)[1]*(x)+(a)[0])
+#define POLYNOM2(x, a) (POLYNOM1((x),(a)+1)*(x)+(a)[0])
+#define POLYNOM3(x, a) (POLYNOM2((x),(a)+1)*(x)+(a)[0])
+#define POLYNOM4(x, a) (POLYNOM3((x),(a)+1)*(x)+(a)[0])
+#define POLYNOM5(x, a) (POLYNOM4((x),(a)+1)*(x)+(a)[0])
+#define POLYNOM6(x, a) (POLYNOM5((x),(a)+1)*(x)+(a)[0])
+#define POLYNOM7(x, a) (POLYNOM6((x),(a)+1)*(x)+(a)[0])
+#define POLYNOM8(x, a) (POLYNOM7((x),(a)+1)*(x)+(a)[0])
+#define POLYNOM9(x, a) (POLYNOM8((x),(a)+1)*(x)+(a)[0])
+#define POLYNOM10(x, a) (POLYNOM9((x),(a)+1)*(x)+(a)[0])
+#define POLYNOM11(x, a) (POLYNOM10((x),(a)+1)*(x)+(a)[0])
+#define POLYNOM12(x, a) (POLYNOM11((x),(a)+1)*(x)+(a)[0])
+#define POLYNOM13(x, a) (POLYNOM12((x),(a)+1)*(x)+(a)[0])
+
+#define M_LN_MAX_D (M_LN2 * DBL_MAX_EXP)
+#define M_LN_MIN_D (M_LN2 * (DBL_MAX_EXP - 1))
--- /dev/null
+/*
+ * (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 <errno.h>
+#include <math.h>
+#include "localmath.h"
+
+
+double
+log(double x)
+{
+ /* log(x) = z*P(z*z)/Q(z*z), z = (x-1)/(x+1), x in [1/sqrt(2), sqrt(2)]
+ */
+ /* 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
+ };
+
+ double z, zsqr;
+ int exponent;
+
+ if (x <= 0) {
+ errno = EDOM;
+ return 0;
+ }
+
+ x = frexp(x, &exponent);
+ while (x < M_1_SQRT2) {
+ x += x;
+ exponent--;
+ }
+ z = (x-1)/(x+1);
+ zsqr = z*z;
+ return z * POLYNOM4(zsqr, p) / POLYNOM4(zsqr, q) + exponent * M_LN2;
+}
--- /dev/null
+/*
+ * (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 <errno.h>
+#include <math.h>
+#include "localmath.h"
+
+double
+log10(double x)
+{
+ if (x <= 0) {
+ errno = EDOM;
+ return 0;
+ }
+
+ return log(x) / M_LN10;
+}
--- /dev/null
+#
+/*
+ * (c) copyright 1987 by the Vrije Universiteit, Amsterdam, The Netherlands.
+ * See the copyright notice in the ACK home directory, in the file "Copyright".
+ */
+/* $Header$ */
+
+ mes 2,EM_WSIZE,EM_PSIZE
+#ifndef NOFLOAT
+ exp $modf
+ pro $modf,0
+ lal 0
+ loi EM_DSIZE
+ loc 1
+ loc EM_WSIZE
+ loc EM_DSIZE
+ cif
+ fif EM_DSIZE
+ lal EM_DSIZE
+ loi EM_PSIZE
+ sti EM_DSIZE
+ ret EM_DSIZE
+ end
+#endif
--- /dev/null
+/*
+ * (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 <errno.h>
+#include <math.h>
+
+double
+pow(double x, double y)
+{
+ double dummy;
+
+ if ((x == 0 && y == 0) ||
+ (x < 0 && modf(y, &dummy) != 0)) {
+ errno = EDOM;
+ return 0;
+ }
+
+ 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;
+ }
+ return val;
+ }
+
+ return exp(log(x) * y);
+}
--- /dev/null
+/*
+ * (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 <errno.h>
+#include <math.h>
+#include "localmath.h"
+
+static double
+sinus(double x, int quadrant)
+{
+ /* 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
+ };
+
+ 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;
+ }
+ 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 (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 {
+ 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) {
+ x = -x;
+ }
+ xsqr = x * x;
+ x = x * POLYNOM5(xsqr, p) / POLYNOM5(xsqr, q);
+ return x;
+}
+
+double
+sin(double x)
+{
+ return sinus(x, 0);
+}
+
+double
+cos(double x)
+{
+ if (x < 0) x = -x;
+ return sinus(x, 1);
+}
--- /dev/null
+/*
+ * (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 <errno.h>
+#include <math.h>
+#include "localmath.h"
+
+double
+sinh(double x)
+{
+ int negx = x < 0;
+
+ if (negx) {
+ x = -x;
+ }
+ if (x > M_LN_MAX_D) {
+ /* exp(x) would overflow */
+ if (x >= M_LN_MAX_D + M_LN2) {
+ /* not representable */
+ x = HUGE_VAL;
+ errno = ERANGE;
+ }
+ else x = exp (x - M_LN2);
+ }
+ else {
+ double expx = exp(x);
+ x = 0.5 * (expx - 1.0/expx);
+ }
+ if (negx) {
+ return -x;
+ }
+ return x;
+}
--- /dev/null
+/*
+ * (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 <errno.h>
+#include <math.h>
+
+#define NITER 5
+
+double
+sqrt(double x)
+{
+ int exponent;
+ double val;
+
+ if (x <= 0) {
+ if (x < 0) errno = EDOM;
+ return 0;
+ }
+
+ val = frexp(x, &exponent);
+ if (exponent & 1) {
+ exponent--;
+ val *= 2;
+ }
+ 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;
+ }
+ return val;
+}
--- /dev/null
+/*
+ * (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 <errno.h>
+#include <math.h>
+#include "localmath.h"
+
+
+double
+tan(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)
+ */
+
+ 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;
+
+ 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;
+ x += x2;
+ x -= n * A2;
+#undef A1
+#undef A2
+ }
+ else {
+ 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_VAL;
+ }
+ else tmp = tmp2 / tmp1;
+ }
+
+ return negative ? -tmp : tmp;
+}
--- /dev/null
+/*
+ * (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 <errno.h>
+#include <math.h>
+#include "localmath.h"
+
+double
+tanh(double x)
+{
+ if (x <= 0.5*M_LN_MIN_D) {
+ return -1;
+ }
+ if (x >= 0.5*M_LN_MAX_D) {
+ return 1;
+ }
+ x = exp(x + x);
+ return (x - 1.0)/(x + 1.0);
+}
public / ack.git / commitdiff