From a18fcb9048f1d777fdb3858d68111b5bd6ed32cf Mon Sep 17 00:00:00 2001 From: ceriel Date: Fri, 22 Jul 1988 16:53:29 +0000 Subject: [PATCH] Initial revision --- lang/cem/libcc/math/LIST | 22 ++++ lang/cem/libcc/math/asin.c | 53 ++++++++++ lang/cem/libcc/math/atan.c | 103 ++++++++++++++++++ lang/cem/libcc/math/atan2.c | 46 ++++++++ lang/cem/libcc/math/ceil.c | 21 ++++ lang/cem/libcc/math/cosh.c | 38 +++++++ lang/cem/libcc/math/exp.c | 67 ++++++++++++ lang/cem/libcc/math/fabs.c | 15 +++ lang/cem/libcc/math/floor.c | 21 ++++ lang/cem/libcc/math/gamma.c | 137 ++++++++++++++++++++++++ lang/cem/libcc/math/hypot.c | 39 +++++++ lang/cem/libcc/math/j0.c | 203 +++++++++++++++++++++++++++++++++++ lang/cem/libcc/math/j1.c | 206 ++++++++++++++++++++++++++++++++++++ lang/cem/libcc/math/jn.c | 121 +++++++++++++++++++++ lang/cem/libcc/math/log.c | 56 ++++++++++ lang/cem/libcc/math/log10.c | 27 +++++ lang/cem/libcc/math/pow.c | 40 +++++++ lang/cem/libcc/math/sin.c | 115 ++++++++++++++++++++ lang/cem/libcc/math/sinh.c | 42 ++++++++ lang/cem/libcc/math/sqrt.c | 41 +++++++ lang/cem/libcc/math/tan.c | 126 ++++++++++++++++++++++ lang/cem/libcc/math/tanh.c | 27 +++++ lang/cem/libcc/math/test.c | 193 +++++++++++++++++++++++++++++++++ 23 files changed, 1759 insertions(+) create mode 100644 lang/cem/libcc/math/LIST create mode 100644 lang/cem/libcc/math/asin.c create mode 100644 lang/cem/libcc/math/atan.c create mode 100644 lang/cem/libcc/math/atan2.c create mode 100644 lang/cem/libcc/math/ceil.c create mode 100644 lang/cem/libcc/math/cosh.c create mode 100644 lang/cem/libcc/math/exp.c create mode 100644 lang/cem/libcc/math/fabs.c create mode 100644 lang/cem/libcc/math/floor.c create mode 100644 lang/cem/libcc/math/gamma.c create mode 100644 lang/cem/libcc/math/hypot.c create mode 100644 lang/cem/libcc/math/j0.c create mode 100644 lang/cem/libcc/math/j1.c create mode 100644 lang/cem/libcc/math/jn.c create mode 100644 lang/cem/libcc/math/log.c create mode 100644 lang/cem/libcc/math/log10.c create mode 100644 lang/cem/libcc/math/pow.c create mode 100644 lang/cem/libcc/math/sin.c create mode 100644 lang/cem/libcc/math/sinh.c create mode 100644 lang/cem/libcc/math/sqrt.c create mode 100644 lang/cem/libcc/math/tan.c create mode 100644 lang/cem/libcc/math/tanh.c create mode 100644 lang/cem/libcc/math/test.c diff --git a/lang/cem/libcc/math/LIST b/lang/cem/libcc/math/LIST new file mode 100644 index 000000000..7e0b35ab8 --- /dev/null +++ b/lang/cem/libcc/math/LIST @@ -0,0 +1,22 @@ +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 diff --git a/lang/cem/libcc/math/asin.c b/lang/cem/libcc/math/asin.c new file mode 100644 index 000000000..0e4b3f05f --- /dev/null +++ b/lang/cem/libcc/math/asin.c @@ -0,0 +1,53 @@ +/* + * (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 + +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); +} diff --git a/lang/cem/libcc/math/atan.c b/lang/cem/libcc/math/atan.c new file mode 100644 index 000000000..787d9ffbe --- /dev/null +++ b/lang/cem/libcc/math/atan.c @@ -0,0 +1,103 @@ +/* + * (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 + +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; +} diff --git a/lang/cem/libcc/math/atan2.c b/lang/cem/libcc/math/atan2.c new file mode 100644 index 000000000..83e6d4aa8 --- /dev/null +++ b/lang/cem/libcc/math/atan2.c @@ -0,0 +1,46 @@ +/* + * (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 + +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; +} diff --git a/lang/cem/libcc/math/ceil.c b/lang/cem/libcc/math/ceil.c new file mode 100644 index 000000000..7d7bb7ea4 --- /dev/null +++ b/lang/cem/libcc/math/ceil.c @@ -0,0 +1,21 @@ +/* + * (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 + */ +} diff --git a/lang/cem/libcc/math/cosh.c b/lang/cem/libcc/math/cosh.c new file mode 100644 index 000000000..4cecd4e5f --- /dev/null +++ b/lang/cem/libcc/math/cosh.c @@ -0,0 +1,38 @@ +/* + * (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 + +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; +} diff --git a/lang/cem/libcc/math/exp.c b/lang/cem/libcc/math/exp.c new file mode 100644 index 000000000..737ae234c --- /dev/null +++ b/lang/cem/libcc/math/exp.c @@ -0,0 +1,67 @@ +/* + * (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 + +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; +} diff --git a/lang/cem/libcc/math/fabs.c b/lang/cem/libcc/math/fabs.c new file mode 100644 index 000000000..26e08336d --- /dev/null +++ b/lang/cem/libcc/math/fabs.c @@ -0,0 +1,15 @@ +/* + * (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; +} diff --git a/lang/cem/libcc/math/floor.c b/lang/cem/libcc/math/floor.c new file mode 100644 index 000000000..e6b9ac652 --- /dev/null +++ b/lang/cem/libcc/math/floor.c @@ -0,0 +1,21 @@ +/* + * (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 + */ +} diff --git a/lang/cem/libcc/math/gamma.c b/lang/cem/libcc/math/gamma.c new file mode 100644 index 000000000..577a6a34f --- /dev/null +++ b/lang/cem/libcc/math/gamma.c @@ -0,0 +1,137 @@ +/* + * (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 + +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)); +} diff --git a/lang/cem/libcc/math/hypot.c b/lang/cem/libcc/math/hypot.c new file mode 100644 index 000000000..1fa6e30aa --- /dev/null +++ b/lang/cem/libcc/math/hypot.c @@ -0,0 +1,39 @@ +/* + * (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); +} diff --git a/lang/cem/libcc/math/j0.c b/lang/cem/libcc/math/j0.c new file mode 100644 index 000000000..3cf93c686 --- /dev/null +++ b/lang/cem/libcc/math/j0.c @@ -0,0 +1,203 @@ +/* + * (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 + + +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); +} diff --git a/lang/cem/libcc/math/j1.c b/lang/cem/libcc/math/j1.c new file mode 100644 index 000000000..71fab5a3a --- /dev/null +++ b/lang/cem/libcc/math/j1.c @@ -0,0 +1,206 @@ +/* + * (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 + + +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); +} diff --git a/lang/cem/libcc/math/jn.c b/lang/cem/libcc/math/jn.c new file mode 100644 index 000000000..d1c60f423 --- /dev/null +++ b/lang/cem/libcc/math/jn.c @@ -0,0 +1,121 @@ +/* + * (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 + +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; + } +} diff --git a/lang/cem/libcc/math/log.c b/lang/cem/libcc/math/log.c new file mode 100644 index 000000000..1b2dc37a2 --- /dev/null +++ b/lang/cem/libcc/math/log.c @@ -0,0 +1,56 @@ +/* + * (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 + +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; +} diff --git a/lang/cem/libcc/math/log10.c b/lang/cem/libcc/math/log10.c new file mode 100644 index 000000000..b2e7cacd7 --- /dev/null +++ b/lang/cem/libcc/math/log10.c @@ -0,0 +1,27 @@ +/* + * (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 + +extern int errno; + +double +log10(x) + double x; +{ + extern double log(); + + if (x <= 0) { + errno = EDOM; + return 0; + } + + return log(x) / M_LN10; +} diff --git a/lang/cem/libcc/math/pow.c b/lang/cem/libcc/math/pow.c new file mode 100644 index 000000000..893f82c4b --- /dev/null +++ b/lang/cem/libcc/math/pow.c @@ -0,0 +1,40 @@ +/* + * (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 + +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); +} diff --git a/lang/cem/libcc/math/sin.c b/lang/cem/libcc/math/sin.c new file mode 100644 index 000000000..97fbdafb4 --- /dev/null +++ b/lang/cem/libcc/math/sin.c @@ -0,0 +1,115 @@ +/* + * (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 + +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); +} diff --git a/lang/cem/libcc/math/sinh.c b/lang/cem/libcc/math/sinh.c new file mode 100644 index 000000000..f5f94b745 --- /dev/null +++ b/lang/cem/libcc/math/sinh.c @@ -0,0 +1,42 @@ +/* + * (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 + +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; +} diff --git a/lang/cem/libcc/math/sqrt.c b/lang/cem/libcc/math/sqrt.c new file mode 100644 index 000000000..4369eaeb4 --- /dev/null +++ b/lang/cem/libcc/math/sqrt.c @@ -0,0 +1,41 @@ +/* + * (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 + +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; +} diff --git a/lang/cem/libcc/math/tan.c b/lang/cem/libcc/math/tan.c new file mode 100644 index 000000000..350cea46d --- /dev/null +++ b/lang/cem/libcc/math/tan.c @@ -0,0 +1,126 @@ +/* + * (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 + +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; +} diff --git a/lang/cem/libcc/math/tanh.c b/lang/cem/libcc/math/tanh.c new file mode 100644 index 000000000..3bceb37e9 --- /dev/null +++ b/lang/cem/libcc/math/tanh.c @@ -0,0 +1,27 @@ +/* + * (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 + +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); +} diff --git a/lang/cem/libcc/math/test.c b/lang/cem/libcc/math/test.c new file mode 100644 index 000000000..a12f19f33 --- /dev/null +++ b/lang/cem/libcc/math/test.c @@ -0,0 +1,193 @@ +/* + * (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 +#include + +#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); + } + } +} -- 2.34.1