From d2f7f252b2be99a7b0f3b185b6cb54ec71bd692e Mon Sep 17 00:00:00 2001 From: eck Date: Wed, 10 May 1989 16:08:14 +0000 Subject: [PATCH] Initial revision --- lang/cem/libcc.ansi/math/asin.c | 47 +++++++++++ lang/cem/libcc.ansi/math/atan.c | 104 +++++++++++++++++++++++ lang/cem/libcc.ansi/math/atan2.c | 42 +++++++++ lang/cem/libcc.ansi/math/ceil.c | 20 +++++ lang/cem/libcc.ansi/math/cosh.c | 33 ++++++++ lang/cem/libcc.ansi/math/exp.c | 65 ++++++++++++++ lang/cem/libcc.ansi/math/fabs.c | 13 +++ lang/cem/libcc.ansi/math/floor.c | 20 +++++ lang/cem/libcc.ansi/math/frexp.e | 20 +++++ lang/cem/libcc.ansi/math/ldexp.c | 36 ++++++++ lang/cem/libcc.ansi/math/localmath.h | 42 +++++++++ lang/cem/libcc.ansi/math/log.c | 53 ++++++++++++ lang/cem/libcc.ansi/math/log10.c | 22 +++++ lang/cem/libcc.ansi/math/modf.e | 24 ++++++ lang/cem/libcc.ansi/math/pow.c | 35 ++++++++ lang/cem/libcc.ansi/math/sin.c | 109 ++++++++++++++++++++++++ lang/cem/libcc.ansi/math/sinh.c | 38 +++++++++ lang/cem/libcc.ansi/math/sqrt.c | 36 ++++++++ lang/cem/libcc.ansi/math/tan.c | 122 +++++++++++++++++++++++++++ lang/cem/libcc.ansi/math/tanh.c | 24 ++++++ 20 files changed, 905 insertions(+) create mode 100644 lang/cem/libcc.ansi/math/asin.c create mode 100644 lang/cem/libcc.ansi/math/atan.c create mode 100644 lang/cem/libcc.ansi/math/atan2.c create mode 100644 lang/cem/libcc.ansi/math/ceil.c create mode 100644 lang/cem/libcc.ansi/math/cosh.c create mode 100644 lang/cem/libcc.ansi/math/exp.c create mode 100644 lang/cem/libcc.ansi/math/fabs.c create mode 100644 lang/cem/libcc.ansi/math/floor.c create mode 100644 lang/cem/libcc.ansi/math/frexp.e create mode 100644 lang/cem/libcc.ansi/math/ldexp.c create mode 100644 lang/cem/libcc.ansi/math/localmath.h create mode 100644 lang/cem/libcc.ansi/math/log.c create mode 100644 lang/cem/libcc.ansi/math/log10.c create mode 100644 lang/cem/libcc.ansi/math/modf.e create mode 100644 lang/cem/libcc.ansi/math/pow.c create mode 100644 lang/cem/libcc.ansi/math/sin.c create mode 100644 lang/cem/libcc.ansi/math/sinh.c create mode 100644 lang/cem/libcc.ansi/math/sqrt.c create mode 100644 lang/cem/libcc.ansi/math/tan.c create mode 100644 lang/cem/libcc.ansi/math/tanh.c diff --git a/lang/cem/libcc.ansi/math/asin.c b/lang/cem/libcc.ansi/math/asin.c new file mode 100644 index 000000000..f5e5a15e5 --- /dev/null +++ b/lang/cem/libcc.ansi/math/asin.c @@ -0,0 +1,47 @@ +/* + * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + * See the copyright notice in the ACK home directory, in the file "Copyright". + * + * Author: Ceriel J.H. Jacobs + */ +/* $Header$ */ + +#include +#include +#include "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); +} diff --git a/lang/cem/libcc.ansi/math/atan.c b/lang/cem/libcc.ansi/math/atan.c new file mode 100644 index 000000000..80de677f7 --- /dev/null +++ b/lang/cem/libcc.ansi/math/atan.c @@ -0,0 +1,104 @@ +/* + * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + * See the copyright notice in the ACK home directory, in the file "Copyright". + * + * Author: Ceriel J.H. Jacobs + */ +/* $Header$ */ + +#include +#include +#include +#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; +} + diff --git a/lang/cem/libcc.ansi/math/atan2.c b/lang/cem/libcc.ansi/math/atan2.c new file mode 100644 index 000000000..541e91746 --- /dev/null +++ b/lang/cem/libcc.ansi/math/atan2.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 +#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; +} diff --git a/lang/cem/libcc.ansi/math/ceil.c b/lang/cem/libcc.ansi/math/ceil.c new file mode 100644 index 000000000..15f649288 --- /dev/null +++ b/lang/cem/libcc.ansi/math/ceil.c @@ -0,0 +1,20 @@ +/* + * (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 + +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 + */ +} diff --git a/lang/cem/libcc.ansi/math/cosh.c b/lang/cem/libcc.ansi/math/cosh.c new file mode 100644 index 000000000..89456e765 --- /dev/null +++ b/lang/cem/libcc.ansi/math/cosh.c @@ -0,0 +1,33 @@ +/* + * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + * See the copyright notice in the ACK home directory, in the file "Copyright". + * + * Author: Ceriel J.H. Jacobs + */ +/* $Header$ */ + +#include +#include +#include "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; +} diff --git a/lang/cem/libcc.ansi/math/exp.c b/lang/cem/libcc.ansi/math/exp.c new file mode 100644 index 000000000..d731d6540 --- /dev/null +++ b/lang/cem/libcc.ansi/math/exp.c @@ -0,0 +1,65 @@ +/* + * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + * See the copyright notice in the ACK home directory, in the file "Copyright". + * + * Author: Ceriel J.H. Jacobs + */ +/* $Header$ */ + +#include +#include +#include +#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; +} diff --git a/lang/cem/libcc.ansi/math/fabs.c b/lang/cem/libcc.ansi/math/fabs.c new file mode 100644 index 000000000..8c190c2c6 --- /dev/null +++ b/lang/cem/libcc.ansi/math/fabs.c @@ -0,0 +1,13 @@ +/* + * (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; +} diff --git a/lang/cem/libcc.ansi/math/floor.c b/lang/cem/libcc.ansi/math/floor.c new file mode 100644 index 000000000..ee0e6de05 --- /dev/null +++ b/lang/cem/libcc.ansi/math/floor.c @@ -0,0 +1,20 @@ +/* + * (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 + +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 + */ +} diff --git a/lang/cem/libcc.ansi/math/frexp.e b/lang/cem/libcc.ansi/math/frexp.e new file mode 100644 index 000000000..922470c60 --- /dev/null +++ b/lang/cem/libcc.ansi/math/frexp.e @@ -0,0 +1,20 @@ +# +/* + * (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 diff --git a/lang/cem/libcc.ansi/math/ldexp.c b/lang/cem/libcc.ansi/math/ldexp.c new file mode 100644 index 000000000..4a48f1d5b --- /dev/null +++ b/lang/cem/libcc.ansi/math/ldexp.c @@ -0,0 +1,36 @@ +/* + * (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 + +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; +} diff --git a/lang/cem/libcc.ansi/math/localmath.h b/lang/cem/libcc.ansi/math/localmath.h new file mode 100644 index 000000000..c73ab5742 --- /dev/null +++ b/lang/cem/libcc.ansi/math/localmath.h @@ -0,0 +1,42 @@ +/* + * 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)) diff --git a/lang/cem/libcc.ansi/math/log.c b/lang/cem/libcc.ansi/math/log.c new file mode 100644 index 000000000..53e709345 --- /dev/null +++ b/lang/cem/libcc.ansi/math/log.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 +#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; +} diff --git a/lang/cem/libcc.ansi/math/log10.c b/lang/cem/libcc.ansi/math/log10.c new file mode 100644 index 000000000..09ab6b882 --- /dev/null +++ b/lang/cem/libcc.ansi/math/log10.c @@ -0,0 +1,22 @@ +/* + * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + * See the copyright notice in the ACK home directory, in the file "Copyright". + * + * Author: Ceriel J.H. Jacobs + */ +/* $Header$ */ + +#include +#include +#include "localmath.h" + +double +log10(double x) +{ + if (x <= 0) { + errno = EDOM; + return 0; + } + + return log(x) / M_LN10; +} diff --git a/lang/cem/libcc.ansi/math/modf.e b/lang/cem/libcc.ansi/math/modf.e new file mode 100644 index 000000000..8396117ba --- /dev/null +++ b/lang/cem/libcc.ansi/math/modf.e @@ -0,0 +1,24 @@ +# +/* + * (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 diff --git a/lang/cem/libcc.ansi/math/pow.c b/lang/cem/libcc.ansi/math/pow.c new file mode 100644 index 000000000..af04e3ea0 --- /dev/null +++ b/lang/cem/libcc.ansi/math/pow.c @@ -0,0 +1,35 @@ +/* + * (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 +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); +} diff --git a/lang/cem/libcc.ansi/math/sin.c b/lang/cem/libcc.ansi/math/sin.c new file mode 100644 index 000000000..de60f99a1 --- /dev/null +++ b/lang/cem/libcc.ansi/math/sin.c @@ -0,0 +1,109 @@ +/* + * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + * See the copyright notice in the ACK home directory, in the file "Copyright". + * + * Author: Ceriel J.H. Jacobs + */ +/* $Header$ */ + +#include +#include +#include "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); +} diff --git a/lang/cem/libcc.ansi/math/sinh.c b/lang/cem/libcc.ansi/math/sinh.c new file mode 100644 index 000000000..39633dc52 --- /dev/null +++ b/lang/cem/libcc.ansi/math/sinh.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 +#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; +} diff --git a/lang/cem/libcc.ansi/math/sqrt.c b/lang/cem/libcc.ansi/math/sqrt.c new file mode 100644 index 000000000..afaf2a490 --- /dev/null +++ b/lang/cem/libcc.ansi/math/sqrt.c @@ -0,0 +1,36 @@ +/* + * (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 + +#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; +} diff --git a/lang/cem/libcc.ansi/math/tan.c b/lang/cem/libcc.ansi/math/tan.c new file mode 100644 index 000000000..3d6911ff1 --- /dev/null +++ b/lang/cem/libcc.ansi/math/tan.c @@ -0,0 +1,122 @@ +/* + * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + * See the copyright notice in the ACK home directory, in the file "Copyright". + * + * Author: Ceriel J.H. Jacobs + */ +/* $Header$ */ + +#include +#include +#include "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; +} diff --git a/lang/cem/libcc.ansi/math/tanh.c b/lang/cem/libcc.ansi/math/tanh.c new file mode 100644 index 000000000..cd1216a38 --- /dev/null +++ b/lang/cem/libcc.ansi/math/tanh.c @@ -0,0 +1,24 @@ +/* + * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. + * See the copyright notice in the ACK home directory, in the file "Copyright". + * + * Author: Ceriel J.H. Jacobs + */ +/* $Header$ */ + +#include +#include +#include "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); +} -- 2.34.1