--- /dev/null
+tail_m.a
+asin.c
+atan2.c
+atan.c
+ceil.c
+cosh.c
+fabs.c
+gamma.c
+hypot.c
+jn.c
+j0.c
+j1.c
+log10.c
+pow.c
+log.c
+sin.c
+sinh.c
+sqrt.c
+tan.c
+tanh.c
+exp.c
+floor.c
--- /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>
+#include <errno.h>
+
+extern int errno;
+
+
+static double
+asin_acos(x, cosfl)
+ double x;
+{
+ int negative = x < 0;
+ extern double sqrt(), atan();
+
+ 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(x)
+ double x;
+{
+ return asin_acos(x, 0);
+}
+
+double
+acos(x)
+ 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 <math.h>
+#include <errno.h>
+
+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)
+ */
+ 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 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 <math.h>
+#include <errno.h>
+
+extern int errno;
+
+double
+atan2(y, x)
+ double x, y;
+{
+ extern double atan();
+ 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$ */
+
+double
+ceil(x)
+ double x;
+{
+ extern double modf();
+ 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 <math.h>
+#include <errno.h>
+
+extern int errno;
+
+double
+cosh(x)
+ double x;
+{
+ extern double exp();
+
+ 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;
+ 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 <math.h>
+#include <errno.h>
+
+extern int errno;
+
+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 */
+
+ 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;
+ extern double floor(), ldexp();
+
+ if (x <= M_LN_MIN_D) {
+ if (x < M_LN_MIN_D) errno = ERANGE;
+ return M_MIN_D;
+ }
+ if (x >= M_LN_MAX_D) {
+ if (x < M_LN_MAX_D) errno = ERANGE;
+ 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;
+ }
+ 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(x)
+ 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$ */
+
+double
+floor(x)
+ double x;
+{
+ extern double modf();
+ 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 <math.h>
+#include <errno.h>
+
+static double
+smallpos_gamma(x)
+ double x;
+{
+ /* Approximation of gamma function using
+ gamma(x) = P(x-2) / Q(x-2) for x in [2,3]
+ */
+ /* Hart & Cheney # 5251 */
+
+ static double p[11] = {
+ -0.2983543278574342138830437659e+06,
+ -0.2384953970018198872468734423e+06,
+ -0.1170494760121780688403854445e+06,
+ -0.3949445048301571936421824091e+05,
+ -0.1046699423827521405330650531e+05,
+ -0.2188218110071816359394795998e+04,
+ -0.3805112208641734657584922631e+03,
+ -0.5283123755635845383718978382e+02,
+ -0.6128571763704498306889428212e+01,
+ -0.5028018054416812467364198750e+00,
+ -0.3343060322330595274515660112e-01
+ };
+
+ static double q[9] = {
+ -0.2983543278574342138830438524e+06,
+ -0.1123558608748644911342306408e+06,
+ 0.5332716689118142157485686311e+05,
+ 0.8571160498907043851961147763e+04,
+ -0.4734865977028211706556819770e+04,
+ 0.1960497612885585838997039621e+03,
+ 0.1257733367869888645966647426e+03,
+ -0.2053126153100672764513929067e+02,
+ 0.1000000000000000000000000000e+01
+ };
+
+ double result = 1.0;
+
+ while (x > 3) {
+ x -= 1.0;
+ result *= x;
+ }
+ while (x < 2) {
+ result /= x;
+ x += 1.0;
+ }
+
+ x -= 2.0;
+
+ return result * POLYNOM10(x, p) / POLYNOM8(x, q);
+}
+
+#define log_sqrt_2pi 0.91893853320467274178032973640561763
+
+int signgam;
+
+static double
+bigpos_loggamma(x)
+ double x;
+{
+ /* computes the log(gamma(x)) function for big arguments
+ using the Stirling form
+ log(gamma(x)) = (x - 0.5)log(x) - x + log(sqrt(2*pi)) + fi(x)
+ where fi(x) = (1/x)*P(1/(x*x))/Q(1/(x*x)) for x in [12,1000]
+ */
+ /* Hart & Cheney # 5468 */
+
+ static double p[4] = {
+ 0.12398282342474941538685913e+00,
+ 0.67082783834332134961461700e+00,
+ 0.64507302912892202513890000e+00,
+ 0.66662907040200752600000000e-01
+ };
+
+ static double q[4] = {
+ 0.14877938810969929846815600e+01,
+ 0.80995271894897557472821400e+01,
+ 0.79966911236636441947720000e+01,
+ 0.10000000000000000000000000e+01
+ };
+
+ double rsq = 1.0/(x*x);
+ extern double log();
+
+ return (x-0.5)*log(x)-x+log_sqrt_2pi+POLYNOM3(rsq, p)/(x*POLYNOM3(rsq, q));
+}
+
+static double
+neg_loggamma(x)
+ double x;
+{
+ /* compute the log(gamma(x)) function for negative values of x,
+ using the rule:
+ -x*gamma(x)*gamma(-x) = pi/sin(z*pi)
+ */
+ extern double sin(), log();
+ double sinpix;
+
+ x = -x;
+ sinpix = sin(M_PI * x);
+ if (sinpix == 0.0) {
+ errno = EDOM;
+ return HUGE;
+ }
+ if (sinpix < 0) sinpix = -sinpix;
+ else signgam = -1;
+ return log(M_PI/(x * smallpos_gamma(x) * sinpix));
+}
+
+double
+gamma(x)
+ double x;
+{
+ /* Wrong name; Actually computes log(gamma(x))
+ */
+ extern double log();
+
+ signgam = 1;
+ if (x <= 0) {
+ return neg_loggamma(x);
+ }
+ if (x > 12.0) {
+ return bigpos_loggamma(x);
+ }
+ return log(smallpos_gamma(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
+hypot(x,y)
+ double x,y;
+{
+ /* Computes sqrt(x*x+y*y), avoiding overflow */
+
+ extern double sqrt();
+
+ if (x < 0) x = -x;
+ if (y < 0) y = -y;
+ if (x > y) {
+ double t = y;
+ y = x;
+ x = t;
+ }
+ /* sqrt(x*x+y*y) = sqrt(y*y*(x*x/(y*y)+1.0)) = y*sqrt(x*x/(y*y)+1.0) */
+ x /= y;
+ return y*sqrt(x*x+1.0);
+}
+
+struct complex {
+ double r,i;
+};
+
+double
+cabs(p_compl)
+ struct complex p_compl;
+{
+ return hypot(p_compl.r, p_compl.i);
+}
--- /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>
+#include <errno.h>
+
+
+static double
+P0(x)
+ double x;
+{
+ /* P0(x) = P(z*z)/Q(z*z) where z = 8/x, with x >= 8 */
+ /* Hart & Cheney # 6554 */
+
+ static double p[9] = {
+ 0.9999999999999999999999995647e+00,
+ 0.5638253933310769952531889297e+01,
+ 0.1124846237418285392887270013e+02,
+ 0.1009280644639441488899111404e+02,
+ 0.4290591487686900980651458361e+01,
+ 0.8374209971661497198619102718e+00,
+ 0.6702347074465611456598882534e-01,
+ 0.1696260729396856143084502774e-02,
+ 0.6463970103128382090713889584e-05
+ };
+
+ static double q[9] = {
+ 0.9999999999999999999999999999e+00,
+ 0.5639352566123269952531467562e+01,
+ 0.1125463057106955935416066535e+02,
+ 0.1010501892629524191262518048e+02,
+ 0.4301396985171094350444425443e+01,
+ 0.8418926086780046799127094223e+00,
+ 0.6784915305473610998681570734e-01,
+ 0.1754416614608056207958880988e-02,
+ 0.7482977995134121064747276923e-05
+ };
+
+ double zsq = 64.0/(x*x);
+
+ return POLYNOM8(zsq, p) / POLYNOM8(zsq, q);
+}
+
+static double
+Q0(x)
+ double x;
+{
+ /* Q0(x) = z*P(z*z)/Q(z*z) where z = 8/x, x >= 8 */
+ /* Hart & Cheney # 6955 */
+ /* Probably typerror in Hart & Cheney; it sais:
+ Q0(x) = x*P(z*z)/Q(z*z)
+ */
+
+ static double p[9] = {
+ -0.1562499999999999999999995808e-01,
+ -0.1111285583113679178917024959e+00,
+ -0.2877685516355036842789761274e+00,
+ -0.3477683453166454475665803194e+00,
+ -0.2093031978191084473537206358e+00,
+ -0.6209520943730206312601003832e-01,
+ -0.8434508346572023650653353729e-02,
+ -0.4414848186188819989871882393e-03,
+ -0.5768946278415631134804064871e-05
+ };
+
+ static double q[10] = {
+ 0.9999999999999999999999999999e+00,
+ 0.7121383005365046745065850254e+01,
+ 0.1848194194302368046679068851e+02,
+ 0.2242327522435983712994071530e+02,
+ 0.1359286169255959339963319677e+02,
+ 0.4089489268101204780080944780e+01,
+ 0.5722140925672174525430730669e+00,
+ 0.3219814230905924725810683346e-01,
+ 0.5299687475496044642364124073e-03,
+ 0.9423249021001925212258428217e-06
+ };
+
+ double zsq = 64.0/(x*x);
+
+ return (8.0/x) * POLYNOM8(zsq, p) / POLYNOM9(zsq, q);
+}
+
+static double
+smallj0(x)
+ double x;
+{
+ /* J0(x) = P(x*x)/Q(x*x) for x in [0,8] */
+ /* Hart & Cheney # 5852 */
+
+ static double p[10] = {
+ 0.1641556014884554385346147435e+25,
+ -0.3943559664767296636012616471e+24,
+ 0.2172018385924539313982287997e+23,
+ -0.4814859952069817648285245941e+21,
+ 0.5345457598841972345381674607e+19,
+ -0.3301538925689637686465426220e+17,
+ 0.1187390681211042949874031474e+15,
+ -0.2479851167896144439689877514e+12,
+ 0.2803148940831953934479400118e+09,
+ -0.1336625500481224741885945416e+06
+ };
+
+ static double q[10] = {
+ 0.1641556014884554385346137617e+25,
+ 0.1603303724440893273539045602e+23,
+ 0.7913043777646405204323616203e+20,
+ 0.2613165313325153278086066185e+18,
+ 0.6429607918826017759289213100e+15,
+ 0.1237672982083407903483177730e+13,
+ 0.1893012093677918995179541438e+10,
+ 0.2263381356781110003609399116e+07,
+ 0.1974019272727281783930443513e+04,
+ 0.1000000000000000000000000000e+01
+ };
+
+ double xsq = x*x;
+
+ return POLYNOM9(xsq, p) / POLYNOM9(xsq, q);
+}
+
+double
+j0(x)
+ double x;
+{
+ /* Use J0(x) = sqrt(2/(pi*x))*(P0(x)*cos(X0)-Q0(x)*sin(X0))
+ where X0 = x - pi/4 for |x| > 8.
+ Use J0(-x) = J0(x).
+ Use direct approximation of smallj0 for |x| <= 8.
+ */
+ extern double sqrt(), sin(), cos();
+
+ if (x < 0) x = -x;
+ if (x > 8.0) {
+ double X0 = x - M_PI_4;
+ return sqrt(M_2_PI/x)*(P0(x)*cos(X0) - Q0(x)*sin(X0));
+ }
+ return smallj0(x);
+}
+
+static double
+smally0_bar(x)
+ double x;
+{
+ /* Y0(x) = Y0BAR(x)+(2/pi)*J0(x)ln(x)
+ Approximation of Y0BAR for 0 <= x <= 8:
+ Y0BAR(x) = P(x*x)/Q(x*x)
+ Hart & Cheney #6250
+ */
+
+ static double p[14] = {
+ -0.2692670958801060448840356941e+14,
+ 0.6467231173109037044444917683e+14,
+ -0.5563036156275660297303897296e+13,
+ 0.1698403391975239335187832821e+12,
+ -0.2606282788256139370857687880e+10,
+ 0.2352841334491277505699488812e+08,
+ -0.1365184412186963659690851354e+06,
+ 0.5371538422626582142170627457e+03,
+ -0.1478903875146718839145348490e+01,
+ 0.2887840299886172125955719069e-02,
+ -0.3977426824263991024666116123e-05,
+ 0.3738169731655229006655176866e-08,
+ -0.2194460874896856106887900645e-11,
+ 0.6208996973821484304384239393e-15
+ };
+
+ static double q[6] = {
+ 0.3648393301278364629844168660e+15,
+ 0.1698390180526960997295118328e+13,
+ 0.3587111679107612117789088586e+10,
+ 0.4337760840406994515845890005e+07,
+ 0.3037977771964348276793136205e+04,
+ 0.1000000000000000000000000000e+01
+ };
+
+ double xsq = x*x;
+
+ return POLYNOM13(xsq, p) / POLYNOM5(xsq, q);
+}
+
+double
+y0(x)
+ double x;
+{
+ extern double sqrt(), sin(), cos(), log();
+
+ if (x <= 0.0) {
+ errno = EDOM;
+ return -HUGE;
+ }
+ if (x > 8.0) {
+ double X0 = x - M_PI_4;
+ return sqrt(M_2_PI/x) * (P0(x)*sin(X0)+Q0(x)*cos(X0));
+ }
+ return smally0_bar(x) + M_2_PI*j0(x)*log(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>
+#include <errno.h>
+
+
+static double
+P1(x)
+ double x;
+{
+ /* P1(x) = P(z*z)/Q(z*z) where z = 8/x, with x >= 8 */
+ /* Hart & Cheney # 6755 */
+
+ static double p[9] = {
+ 0.1000000000000000000000000489e+01,
+ 0.5581663300347182292169450071e+01,
+ 0.1100186625131173123750501118e+02,
+ 0.9727139359130463694593683431e+01,
+ 0.4060011483142278994462590992e+01,
+ 0.7742832212665311906917358099e+00,
+ 0.6021617752811098752098248630e-01,
+ 0.1482350677236405118074646993e-02,
+ 0.6094215148131061431667573909e-05
+ };
+
+ static double q[9] = {
+ 0.9999999999999999999999999999e+00,
+ 0.5579832245659682292169922224e+01,
+ 0.1099168447731617288972771040e+02,
+ 0.9707206835125961446797916892e+01,
+ 0.4042610016540342097334497865e+01,
+ 0.7671965204303836019508430169e+00,
+ 0.5893258668794493100786371406e-01,
+ 0.1393993644981256852404222530e-02,
+ 0.4585597769784750669754696825e-05
+ };
+
+ double zsq = 64.0/(x*x);
+
+ return POLYNOM8(zsq, p) / POLYNOM8(zsq, q);
+}
+
+static double
+Q1(x)
+ double x;
+{
+ /* Q1(x) = z*P(z*z)/Q(z*z) where z = 8/x, x >= 8 */
+ /* Hart & Cheney # 7157 */
+ /* Probably typerror in Hart & Cheney; it sais:
+ Q1(x) = x*P(z*z)/Q(z*z)
+ */
+
+ static double p[9] = {
+ 0.4687499999999999999999995275e-01,
+ 0.3302394516691663879252493748e+00,
+ 0.8456888491208195767613862428e+00,
+ 0.1008551084218946085420665147e+01,
+ 0.5973407972399900690521296181e+00,
+ 0.1737697433393258207540273097e+00,
+ 0.2303862814819568573893610740e-01,
+ 0.1171224207976250587945594946e-02,
+ 0.1486418220337492918307904804e-04
+ };
+
+ static double q[10] = {
+ 0.9999999999999999999999999999e+00,
+ 0.7049380763213049609070823421e+01,
+ 0.1807129960468949760845562209e+02,
+ 0.2159171174362827330505421695e+02,
+ 0.1283239297740546866114600499e+02,
+ 0.3758349275324260869598403931e+01,
+ 0.5055985453754739528620657666e+00,
+ 0.2665604326323907148063400439e-01,
+ 0.3821140353404633025596424652e-03,
+ 0.3206696590241261037875154062e-06
+ };
+
+ double zsq = 64.0/(x*x);
+
+ return (8.0/x) * POLYNOM8(zsq, p) / POLYNOM9(zsq, q);
+}
+
+static double
+smallj1(x)
+ double x;
+{
+ /* J1(x) = x*P(x*x)/Q(x*x) for x in [0,8] */
+ /* Hart & Cheney # 6054 */
+
+ static double p[10] = {
+ 0.1921176307760798128049021316e+25,
+ -0.2226092031387396254771375773e+24,
+ 0.7894463902082476734673226741e+22,
+ -0.1269424373753606065436561036e+21,
+ 0.1092152214043184787101134641e+19,
+ -0.5454629264396819144157448868e+16,
+ 0.1634659487571284628830445048e+14,
+ -0.2909662785381647825756152444e+11,
+ 0.2853433451054763915026471449e+08,
+ -0.1197705712815379389149134705e+05
+ };
+
+ static double q[10] = {
+ 0.3842352615521596256098041912e+25,
+ 0.3507567066272028105798868716e+23,
+ 0.1611334311633414344007062889e+21,
+ 0.4929612313959850319632645381e+18,
+ 0.1117536965288162684489793105e+16,
+ 0.1969278625584719037168592923e+13,
+ 0.2735606122949877990248154504e+10,
+ 0.2940957355049651347475558106e+07,
+ 0.2274736606126590905134610965e+04,
+ 0.1000000000000000000000000000e+01
+ };
+
+ double xsq = x*x;
+
+ return x * POLYNOM9(xsq, p) / POLYNOM9(xsq, q);
+}
+
+double
+j1(x)
+ double x;
+{
+ /* Use J1(x) = sqrt(2/(pi*x))*(P1(x)*cos(X1)-Q1(x)*sin(X1))
+ where X1 = x - 3*pi/4 for |x| > 8.
+ Use J1(-x) = -J1(x).
+ Use direct approximation of smallj1 for |x| <= 8.
+ */
+ extern double sqrt(), sin(), cos();
+ int negative = x < 0.0;
+
+ if (negative) x = -x;
+ if (x > 8.0) {
+ double X1 = x - (M_PI - M_PI_4);
+ x = sqrt(M_2_PI/x)*(P1(x)*cos(X1) - Q1(x)*sin(X1));
+ }
+ else x = smallj1(x);
+ if (negative) return -x;
+ return x;
+}
+
+static double
+smally1_bar(x)
+ double x;
+{
+ /* Y1(x) = Y1BAR(x)+(2/pi)*(J1(x)ln(x) - 1/x)
+ Approximation of Y1BAR for 0 <= x <= 8:
+ Y1BAR(x) = x*P(x*x)/Q(x*x)
+ Hart & Cheney # 6449
+ */
+
+ static double p[10] = {
+ -0.5862655424363443992938931700e+24,
+ 0.1570668341992328458208364904e+24,
+ -0.7351681299005467428400402479e+22,
+ 0.1390658785759080111485190942e+21,
+ -0.1339544201526785345938109179e+19,
+ 0.7290257386242270629526344379e+16,
+ -0.2340575603057015935501295099e+14,
+ 0.4411516199185230690878878903e+11,
+ -0.4542128738770213026987060358e+08,
+ 0.1988612563465350530472715888e+05
+ };
+
+ static double q[10] = {
+ 0.2990279721605116022908679994e+25,
+ 0.2780285010357803058127175655e+23,
+ 0.1302687474507355553192845146e+21,
+ 0.4071330372239164349602952937e+18,
+ 0.9446611865086570116528399283e+15,
+ 0.1707657951197456205887347694e+13,
+ 0.2440358986882941823431612517e+10,
+ 0.2708852767034077697963790196e+07,
+ 0.2174361138333330803617969305e+04,
+ 0.1000000000000000000000000000e+01
+ };
+
+ double xsq = x*x;
+
+ return x * POLYNOM9(xsq, p) / POLYNOM9(xsq, q);
+}
+
+double
+y1(x)
+ double x;
+{
+ extern double sqrt(), sin(), cos(), log();
+
+ if (x <= 0.0) {
+ errno = EDOM;
+ return -HUGE;
+ }
+ if (x > 8.0) {
+ double X1 = x - (M_PI - M_PI_4);
+ return sqrt(M_2_PI/x) * (P1(x)*sin(X1)+Q1(x)*cos(X1));
+ }
+ return smally1_bar(x) + M_2_PI*(j1(x)*log(x) - 1/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>
+#include <errno.h>
+
+double
+yn(n, x)
+ double x;
+{
+ /* Use y0, y1, and the recurrence relation
+ y(n+1,x) = 2*n*y(n,x)/x - y(n-1, x).
+ According to Hart & Cheney, this is stable for all
+ x, n.
+ Also use: y(-n,x) = (-1)^n * y(n, x)
+ */
+
+ int negative = 0;
+ extern double y0(), y1();
+ double yn1, yn2;
+ register int i;
+
+ if (x <= 0) {
+ errno = EDOM;
+ return -HUGE;
+ }
+
+ if (n < 0) {
+ n = -n;
+ negative = (n % 2);
+ }
+
+ if (n == 0) return y0(x);
+ if (n == 1) return y1(x);
+
+ yn2 = y0(x);
+ yn1 = y1(x);
+ for (i = 1; i < n; i++) {
+ double tmp = yn1;
+ yn1 = (i*2)*yn1/x - yn2;
+ yn2 = tmp;
+ }
+ if (negative) return -yn1;
+ return yn1;
+}
+
+double
+jn(n, x)
+ double x;
+{
+ /* Unfortunately, according to Hart & Cheney, the recurrence
+ j(n+1,x) = 2*n*j(n,x)/x - j(n-1,x) is unstable for
+ increasing n, except when x > n.
+ However, j(n,x)/j(n-1,x) = 2/(2*n-x*x/(2*(n+1)-x*x/( ....
+ (a continued fraction).
+ We can use this to determine KJn and KJn-1, where K is a
+ normalization constant not yet known. This enables us
+ to determine KJn-2, ...., KJ1, KJ0. Now we can use the
+ J0 or J1 approximation to determine K.
+ Use: j(-n, x) = (-1)^n * j(n, x)
+ j(n, -x) = (-1)^n * j(n, x)
+ */
+
+ extern double j0(), j1();
+
+ if (n < 0) {
+ n = -n;
+ x = -x;
+ }
+
+ if (n == 0) return j0(x);
+ if (n == 1) return j1(x);
+ if (x > n) {
+ /* in this case, the recurrence relation is stable for
+ increasing n, so we use that.
+ */
+ double jn2 = j0(x), jn1 = j1(x);
+ register int i;
+
+ for (i = 1; i < n; i++) {
+ double tmp = jn1;
+ jn1 = (2*i)*jn1/x - jn2;
+ jn2 = tmp;
+ }
+ return jn1;
+ }
+ {
+ /* we first compute j(n,x)/j(n-1,x) */
+ register int i;
+ double quotient = 0.0;
+ double xsqr = x*x;
+ double jn1, jn2;
+
+ for (i = 20; /* ??? how many do we need ??? */
+ i > 0; i--) {
+ quotient = xsqr/(2*(i+n) - quotient);
+ }
+ quotient = x / (2*n - quotient);
+
+ jn1 = quotient;
+ jn2 = 1.0;
+ for (i = n-1; i > 0; i--) {
+ /* recurrence relation is stable for decreasing n
+ */
+ double tmp = jn2;
+ jn2 = (2*i)*jn2/x - jn1;
+ jn1 = tmp;
+ }
+ /* So, now we have K*Jn = quotient and K*J0 = jn2.
+ Now it is easy; compute real j0, this gives K = jn2/j0,
+ and this then gives Jn = quotient/K = j0 * quotient / jn2.
+ */
+ return j0(x)*quotient/jn2;
+ }
+}
--- /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>
+#include <errno.h>
+
+extern int errno;
+
+double
+log(x)
+ 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
+ };
+
+ extern double frexp();
+ 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 <math.h>
+#include <errno.h>
+
+extern int errno;
+
+double
+log10(x)
+ double x;
+{
+ extern double log();
+
+ if (x <= 0) {
+ errno = EDOM;
+ return 0;
+ }
+
+ return log(x) / M_LN10;
+}
--- /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>
+#include <errno.h>
+
+extern int errno;
+
+double
+pow(x,y)
+ double x,y;
+{
+ double dummy;
+ extern double modf(), exp(), log();
+
+ 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 <math.h>
+#include <errno.h>
+
+extern int errno;
+
+static double
+sinus(x, quadrant)
+ 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
+ };
+
+ 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 {
+ 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) {
+ x = -x;
+ }
+ xsqr = x * x;
+ x = x * POLYNOM5(xsqr, p) / POLYNOM5(xsqr, q);
+ return x;
+}
+
+double
+sin(x)
+ double x;
+{
+ return sinus(x, 0);
+}
+
+double
+cos(x)
+ 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 <math.h>
+#include <errno.h>
+
+extern int errno;
+
+double
+sinh(x)
+ double x;
+{
+ int negx = x < 0;
+ extern double exp();
+
+ 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;
+ 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 <math.h>
+#include <errno.h>
+
+extern int errno;
+
+#define NITER 5
+
+double
+sqrt(x)
+ double x;
+{
+ extern double frexp(), ldexp();
+ 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 <math.h>
+#include <errno.h>
+
+extern int errno;
+
+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)
+ */
+
+ 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 {
+ 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;
+}
--- /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>
+#include <errno.h>
+
+double
+tanh(x)
+ double x;
+{
+ extern double exp();
+
+ 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);
+}
--- /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
+ */
+
+#include <math.h>
+#include <stdio.h>
+
+#define EPS_D 5.0e-14
+main()
+{
+ testsqrt();
+ testtrig();
+ testexplog();
+ testgamma();
+ testbessel();
+}
+
+dotest(s, x, d, v)
+ char *s;
+ double x, d, v;
+{
+ double fabs();
+
+ if (fabs((v - d) / (fabs(v) < EPS_D ? 1.0 : v)) > EPS_D) {
+ printf(s, x);
+ printf(" = %.16e, should be %.16e\n", d, v);
+ }
+}
+
+testsqrt()
+{
+#define SQRT2 M_SQRT2
+#define SQRT10 3.16227766016837933199889354443271853
+
+ double x, val;
+ extern double sqrt();
+
+ dotest("sqrt(%.1f)", 2.0, sqrt(2.0), SQRT2);
+ dotest("sqrt(%.1f)", 10.0, sqrt(10.0), SQRT10);
+
+ for (x = 0.1; x < 0.1e20; x += x) {
+ val = sqrt(x);
+ dotest("sqrt(%.1f)^2", x, val*val, x);
+ }
+}
+
+testtrig()
+{
+#define SINPI_24 0.13052619222005159154840622789548901
+#define SINPI_16 0.19509032201612826784828486847702224
+#define SINPI_12 0.25881904510252076234889883762404832
+#define SINPI_6 0.5
+#define SINPI_4 M_1_SQRT2
+#define SINPI_3 0.86602540378443864676372317075293618
+#define SINPI_2 1.0
+#define SIN0 0.0
+
+ double x;
+ extern double sin(), cos(), tan(), asin(), acos(), atan(), fabs();
+
+ dotest("sin(0)", 0.0, sin(0.0), SIN0);
+ dotest("sin(pi/24)", M_PI/24 , sin(M_PI/24), SINPI_24);
+ dotest("sin(pi/16)", M_PI/16 , sin(M_PI/16), SINPI_16);
+ dotest("sin(pi/12)", M_PI/12 , sin(M_PI/12), SINPI_12);
+ dotest("sin(pi/6)", M_PI/6 , sin(M_PI/6), SINPI_6);
+ dotest("sin(pi/4)", M_PI_4 , sin(M_PI_4), SINPI_4);
+ dotest("sin(pi/3)", M_PI/3 , sin(M_PI/3), SINPI_3);
+ dotest("sin(pi/2)", M_PI_2 , sin(M_PI_2), SINPI_2);
+ dotest("sin(pi)", 0.0, sin(M_PI), SIN0);
+ dotest("sin(3*pi/2)", 0.0, sin(M_PI+M_PI_2), -SINPI_2);
+
+ dotest("sin(-pi/24)", -M_PI/24 , sin(-M_PI/24), -SINPI_24);
+ dotest("sin(-pi/16)", -M_PI/16 , sin(-M_PI/16), -SINPI_16);
+ dotest("sin(-pi/12)", -M_PI/12 , sin(-M_PI/12), -SINPI_12);
+ dotest("sin(-pi/6)", -M_PI/6 , sin(-M_PI/6), -SINPI_6);
+ dotest("sin(-pi/4)", -M_PI_4 , sin(-M_PI_4), -SINPI_4);
+ dotest("sin(-pi/3)", -M_PI/3 , sin(-M_PI/3), -SINPI_3);
+ dotest("sin(-pi/2)", -M_PI_2 , sin(-M_PI_2), -SINPI_2);
+
+ dotest("cos(pi/2)", M_PI_2, cos(M_PI_2), SIN0);
+ dotest("cos(11pi/24)", M_PI/24 , cos(11*M_PI/24), SINPI_24);
+ dotest("cos(7pi/16)", M_PI/16 , cos(7*M_PI/16), SINPI_16);
+ dotest("cos(5pi/12)", M_PI/12 , cos(5*M_PI/12), SINPI_12);
+ dotest("cos(pi/3)", M_PI/6 , cos(M_PI/3), SINPI_6);
+ dotest("cos(pi/4)", M_PI_4 , cos(M_PI_4), SINPI_4);
+ dotest("cos(pi/6)", M_PI/3 , cos(M_PI/6), SINPI_3);
+ dotest("cos(0)", M_PI_2 , cos(0), SINPI_2);
+ dotest("cos(pi)", M_PI , cos(M_PI), -SINPI_2);
+ dotest("cos(3pi/2)", M_PI , cos(M_PI+M_PI_2), SIN0);
+
+ dotest("cos(-pi/2)", M_PI_2, cos(-M_PI_2), SIN0);
+ dotest("cos(-11pi/24)", M_PI/24 , cos(-11*M_PI/24), SINPI_24);
+ dotest("cos(-7pi/16)", M_PI/16 , cos(-7*M_PI/16), SINPI_16);
+ dotest("cos(-5pi/12)", M_PI/12 , cos(-5*M_PI/12), SINPI_12);
+ dotest("cos(-pi/3)", M_PI/6 , cos(-M_PI/3), SINPI_6);
+ dotest("cos(-pi/4)", M_PI_4 , cos(-M_PI_4), SINPI_4);
+ dotest("cos(-pi/6)", M_PI/3 , cos(-M_PI/6), SINPI_3);
+
+ for (x = -10; x <= 10; x += 0.5) {
+ dotest("sin+2*pi-sin(%.2f)", x, sin(x+M_2PI)-sin(x), 0.0);
+ dotest("cos+2*pi-cos(%.2f)", x, cos(x+M_2PI)-cos(x), 0.0);
+ dotest("tan+2*pi-tan(%.2f)", x, tan(x+M_2PI)-tan(x), 0.0);
+ dotest("tan+pi-tan(%.2f)", x, tan(x+M_PI)-tan(x), 0.0);
+ }
+
+ for (x = -1.5; x <= 1.5; x += 0.1) {
+ dotest("asin(sin(%.2f))", x, asin(sin(x)), x);
+ dotest("acos(cos(%.2f))", x, acos(cos(x)), fabs(x));
+ dotest("atan(tan(%.2f))", x, atan(tan(x)), x);
+ }
+}
+
+testexplog()
+{
+#define EXPMIN1 0.36787944117144232159552377016146087 /* exp(-1) */
+#define EXPMIN1_4 0.77880078307140486824517026697832065 /* exp(-1/4) */
+#define EXP0 1.0 /* exp(0) */
+#define EXP1_4 1.28402541668774148407342056806243646 /* exp(1/4) */
+#define EXP1 M_E /* exp(1) */
+#define LN1 0.0 /* log(1) */
+#define LN2 M_LN2 /* log(2) */
+#define LN4 1.38629436111989061883446424291635313 /* log(4) */
+#define LNE 1.0 /* log(e) */
+#define LN10 M_LN10 /* log(10) */
+
+ extern double exp(), log();
+ double x;
+
+ dotest("exp(%.2f)", -1.0, exp(-1.0), EXPMIN1);
+ dotest("exp(%.2f)", -0.25, exp(-0.25), EXPMIN1_4);
+ dotest("exp(%.2f)", 0.0, exp(0.0), EXP0);
+ dotest("exp(%.2f)", 0.25, exp(0.25), EXP1_4);
+ dotest("exp(%.2f)", 1.0, exp(1.0), EXP1);
+
+ dotest("log(%.2f)", 1.0, log(1.0), LN1);
+ dotest("log(%.2f)", 2.0, log(2.0), LN2);
+ dotest("log(%.2f)", 4.0, log(4.0), LN4);
+ dotest("log(%.2f)", 10.0, log(10.0), LN10);
+ dotest("log(e)", M_E, log(M_E), LNE);
+
+ for (x = -30.0; x <= 30.0; x += 0.5) {
+ dotest("log(exp(%.2f))", x, log(exp(x)), x);
+ }
+}
+
+testgamma()
+{
+ double x, xfac;
+ extern double gamma(), exp();
+
+ for (x = 1.0, xfac = 1.0; x < 30.0; x += 1.0) {
+ dotest("exp(gamma(%.2f))", x, exp(gamma(x)), xfac);
+ xfac *= x;
+ }
+}
+
+testbessel()
+{
+#define J0__PI_4 0.85163191370480801270040601506092607 /* j0(pi/4) */
+#define J0__PI_2 0.47200121576823476744766838787250096 /* j0(pi/2) */
+#define J1__PI_4 0.36318783834686733179559374778892472 /* j1(pi/4) */
+#define J1__PI_2 0.56682408890587393771124496346716028 /* j1(pi/2) */
+#define J10__PI_4 0.00000000002369974904082422018721148 /* j10(p1/4) */
+#define J10__PI_2 0.00000002326614794865976450546482206 /* j10(pi/2) */
+
+ extern double j0(), j1(), jn(), yn();
+ register int n;
+ double x;
+ extern char *sprintf();
+ char buf[100];
+
+ dotest("j0(pi/4)", M_PI_4, j0(M_PI_4), J0__PI_4);
+ dotest("j0(pi/2)", M_PI_2, j0(M_PI_2), J0__PI_2);
+ dotest("j1(pi/4)", M_PI_4, j1(M_PI_4), J1__PI_4);
+ dotest("j1(pi/2)", M_PI_2, j1(M_PI_2), J1__PI_2);
+ dotest("j10(pi/4)", M_PI_4, jn(10,M_PI_4), J10__PI_4);
+ dotest("j10(pi/2)", M_PI_2, jn(10,M_PI_2), J10__PI_2);
+
+ /* Also check consistency using the Wronskian relation
+ jn(n+1,x)*yn(n, x) - jn(n,x)*yn(n+1,x) = 2/(pi*x)
+ */
+
+ for (x = 0.1; x < 20.0; x += 0.5) {
+ double two_over_pix = M_2_PI/x;
+
+ for (n = 0; n <= 10; n++) {
+ dotest(sprintf(buf, "jn(%d,%.2f)*yn(%d,%.2f)-jn(%d,%.2f)*yn(%d,%.2f)",n+1,x,n,x,n,x,n+1,x), x, jn(n+1,x)*yn(n,x)-jn(n,x)*yn(n+1,x),M_2_PI/x);
+ }
+ }
+}
public / ack.git / commitdiff