d0p1/ack
1
0
Fork 0
Code Issues Pull requests Packages Projects Releases Wiki Activity

changed from Hart & Cheney to Cody & Waite

Browse source
This commit is contained in:
eck 1989-12-18 15:44:36 +00:00
parent 6e2b44962f
commit d43142d811
19 changed files with 483 additions and 374 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

22
lang/cem/libcc.ansi/math/.distr Normal file
View file

@ -0,0 +1,22 @@
LIST
Makefile
asin.c
atan.c
atan2.c
ceil.c
exp.c
fabs.c
floor.c
fmod.c
frexp.e
ldexp.c
localmath.h
log.c
log10.c
modf.e
pow.c
sin.c
sinh.c
sqrt.c
tan.c
tanh.c

7
lang/cem/libcc.ansi/math/LIST
View file

@ -1,10 +1,8 @@
mlib
localmath.h
asin.c
atan2.c
atan.c
ceil.c
cosh.c
fabs.c
pow.c
log10.c
@ -12,8 +10,11 @@ log.c
sin.c
sinh.c
sqrt.c
ldexp.c
tan.c
tanh.c
exp.c
ldexp.c
fmod.c
floor.c
frexp.e
modf.e

31
lang/cem/libcc.ansi/math/Makefile Normal file
View file

@ -0,0 +1,31 @@
CFLAGS=-L -LIB
.SUFFIXES: .o .e .c
.e.o:
$(CC) $(CFLAGS) -c -o $@ $*.e
clean:
rm -rf asin.o atan2.o atan.o ceil.o fabs.o pow.o log10.o \
log.o sin.o sinh.o sqrt.o tan.o tanh.o exp.o ldexp.o \
fmod.o floor.o frexp.o modf.o OLIST
asin.o:
atan2.o:
atan.o:
ceil.o:
fabs.o:
pow.o:
log10.o:
log.o:
sin.o:
sinh.o:
sqrt.o:
tan.o:
tanh.o:
exp.o:
ldexp.o:
fmod.o:
floor.o:
frexp.o:
modf.o:

50
lang/cem/libcc.ansi/math/asin.c
View file

@ -6,31 +6,61 @@
*/
/* $Header$ */
#include <errno.h>
#include <math.h>
#include <errno.h>
#include "localmath.h"
static double
asin_acos(double x, int cosfl)
{
int negative = x < 0;
int i;
double g;
static double p[] = {
-0.27368494524164255994e+2,
0.57208227877891731407e+2,
-0.39688862997540877339e+2,
0.10152522233806463645e+2,
-0.69674573447350646411e+0
};
static double q[] = {
-0.16421096714498560795e+3,
0.41714430248260412556e+3,
-0.38186303361750149284e+3,
0.15095270841030604719e+3,
-0.23823859153670238830e+2,
1.0
};
if (negative) {
x = -x;
}
if (x > 1) {
errno = EDOM;
return 0;
if (x > 0.5) {
i = 1;
if (x > 1) {
errno = EDOM;
return 0;
}
g = 0.5 - 0.5 * x;
x = - sqrt(g);
x += x;
}
if (x == 1) {
x = M_PI_2;
else {
/* ??? avoid underflow ??? */
i = 0;
g = x * x;
}
else x = atan(x/sqrt(1-x*x));
if (negative) x = -x;
x += x * g * POLYNOM4(g, p) / POLYNOM5(g, q);
if (cosfl) {
return M_PI_2 - x;
if (! negative) x = -x;
}
if ((cosfl == 0) == (i == 1)) {
x = (x + M_PI_4) + M_PI_4;
}
else if (cosfl && negative && i == 1) {
x = (x + M_PI_2) + M_PI_2;
}
if (! cosfl && negative) x = -x;
return x;
}

119
lang/cem/libcc.ansi/math/atan.c
View file

@ -6,7 +6,6 @@
*/
/* $Header$ */
#include <errno.h>
#include <float.h>
#include <math.h>
#include "localmath.h"
@ -14,91 +13,55 @@
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)
/* Algorithm and coefficients from:
"Software manual for the elementary functions"
by W.J. Cody and W. Waite, Prentice-Hall, 1980
*/
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 p[] = {
-0.13688768894191926929e+2,
-0.20505855195861651981e+2,
-0.84946240351320683534e+1,
-0.83758299368150059274e+0
};
static double q[] = {
0.41066306682575781263e+2,
0.86157349597130242515e+2,
0.59578436142597344465e+2,
0.15024001160028576121e+2,
1.0
};
static double a[] = {
0.0,
0.52359877559829887307710723554658381, /* pi/6 */
M_PI_2,
1.04719755119659774615421446109316763 /* pi/3 */
};
static double q[6] = {
0.7698297257888171026986294911e+03,
0.1813892701754635858982709369e+04,
0.1484049607102276827437401170e+04,
0.4904645326203706217748848797e+03,
0.5593479839280348664778328000e+02,
0.1000000000000000000000000000e+01
};
int neg = x < 0;
int n;
double g;
int negative = x < 0.0;
register struct precomputed *pr = prec;
if (negative) {
if (neg) {
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));
if (x > 1.0) {
x = 1.0/x;
n = 2;
}
else {
double xsq = x*x;
else n = 0;
x = x * POLYNOM4(xsq, p)/POLYNOM5(xsq, q);
if (x > 0.26794919243112270647) { /* 2-sqtr(3) */
n = n + 1;
x = (((0.73205080756887729353*x-0.5)-0.5)+x)/
(1.73205080756887729353+x);
}
return negative ? -x : x;
/* ??? avoid underflow ??? */
g = x * x;
x += x * g * POLYNOM3(g, p) / POLYNOM4(g, q);
if (n > 1) x = -x;
x += a[n];
return neg ? -x : x;
}

2
lang/cem/libcc.ansi/math/atan2.c
View file

@ -6,8 +6,8 @@
*/
/* $Header$ */
#include <errno.h>
#include <math.h>
#include <errno.h>
#include "localmath.h"
double

69
lang/cem/libcc.ansi/math/exp.c
View file

@ -6,34 +6,35 @@
*/
/* $Header$ */
#include <errno.h>
#include <float.h>
#include <math.h>
#include <float.h>
#include <errno.h>
#include "localmath.h"
double
exp(double x)
{
/* 2**x = (Q(x*x)+x*P(x*x))/(Q(x*x)-x*P(x*x)) for x in [0,0.5] */
/* Hart & Cheney #1069 */
/* Algorithm and coefficients from:
"Software manual for the elementary functions"
by W.J. Cody and W. Waite, Prentice-Hall, 1980
*/
static double p[3] = {
0.2080384346694663001443843411e+07,
0.3028697169744036299076048876e+05,
0.6061485330061080841615584556e+02
static double p[] = {
0.25000000000000000000e+0,
0.75753180159422776666e-2,
0.31555192765684646356e-4
};
static double q[4] = {
0.6002720360238832528230907598e+07,
0.3277251518082914423057964422e+06,
0.1749287689093076403844945335e+04,
0.1000000000000000000000000000e+01
static double q[] = {
0.50000000000000000000e+0,
0.56817302698551221787e-1,
0.63121894374398503557e-3,
0.75104028399870046114e-6
};
int negative = x < 0;
int ipart, large = 0;
double xsqr, xPxx, Qxx;
double xn, g;
int n;
int negative = x < 0;
if (x <= M_LN_MIN_D) {
if (x < M_LN_MIN_D) errno = ERANGE;
@ -44,22 +45,24 @@ exp(double x)
return DBL_MAX;
}
if (negative) x = -x;
/* ??? avoid underflow ??? */
n = x * M_LOG2E + 0.5; /* 1/ln(2) = log2(e), 0.5 added for rounding */
xn = n;
{
double x1 = (long) x;
double x2 = x - x1;
g = ((x1-xn*0.693359375)+x2) - xn*(-2.1219444005469058277e-4);
}
if (negative) {
x = -x;
g = -g;
n = -n;
}
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;
xn = g * g;
x = g * POLYNOM2(xn, p);
n += 1;
return (ldexp(0.5 + x/(POLYNOM3(xn, q) - x), n));
}

34
lang/cem/libcc.ansi/math/fmod.c Normal file
View file

@ -0,0 +1,34 @@
/*
* (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
* See the copyright notice in the ACK home directory, in the file "Copyright".
*
* Author: Hans van Eck
*/
/* $Header$ */
#include <math.h>
#include <errno.h>
double
fmod(double x, double y)
{
long i;
double val;
double frac;
if (y == 0) {
errno = EDOM;
return 0;
}
frac = modf( x / y, &val);
return frac * y;
/*
val = x / y;
if (val > LONG_MIN && val < LONG_MAX) {
i = val;
return x - i * y;
}
*/
}

14
lang/cem/libcc.ansi/math/frexp.e
View file

@ -5,16 +5,16 @@
*/
/* $Header$ */
mes 2,EM_WSIZE,EM_PSIZE
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
loi _EM_DSIZE
fef _EM_DSIZE
lal _EM_DSIZE
loi _EM_PSIZE
sti _EM_WSIZE
ret _EM_DSIZE
end
#endif

2
lang/cem/libcc.ansi/math/localmath.h
View file

@ -39,4 +39,4 @@
#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))
#define M_LN_MIN_D (M_LN2 * (DBL_MIN_EXP - 1))

62
lang/cem/libcc.ansi/math/log.c
View file

@ -6,48 +6,54 @@
*/
/* $Header$ */
#include <errno.h>
#include <math.h>
#include <errno.h>
#include "localmath.h"
double
log(double x)
{
/* log(x) = z*P(z*z)/Q(z*z), z = (x-1)/(x+1), x in [1/sqrt(2), sqrt(2)]
/* Algorithm and coefficients from:
"Software manual for the elementary functions"
by W.J. Cody and W. Waite, Prentice-Hall, 1980
*/
/* Hart & Cheney #2707 */
static double p[5] = {
0.7504094990777122217455611007e+02,
-0.1345669115050430235318253537e+03,
0.7413719213248602512779336470e+02,
-0.1277249755012330819984385000e+02,
0.3327108381087686938144000000e+00
static double a[] = {
-0.64124943423745581147e2,
0.16383943563021534222e2,
-0.78956112887491257267e0
};
static double b[] = {
-0.76949932108494879777e3,
0.31203222091924532844e3,
-0.35667977739034646171e2,
1.0
};
static double q[5] = {
0.3752047495388561108727775374e+02,
-0.7979028073715004879439951583e+02,
0.5616126132118257292058560360e+02,
-0.1450868091858082685362325000e+02,
0.1000000000000000000000000000e+01
};
double znum, zden, z, w;
int exponent;
double z, zsqr;
int exponent;
if (x <= 0) {
if (x < 0) {
errno = EDOM;
return 0;
return -HUGE_VAL;
}
else if (x == 0) {
errno = ERANGE;
return -HUGE_VAL;
}
x = frexp(x, &exponent);
while (x < M_1_SQRT2) {
x += x;
if (x > M_1_SQRT2) {
znum = (x - 0.5) - 0.5;
zden = x * 0.5 + 0.5;
}
else {
znum = x - 0.5;
zden = znum * 0.5 + 0.5;
exponent--;
}
z = (x-1)/(x+1);
zsqr = z*z;
return z * POLYNOM4(zsqr, p) / POLYNOM4(zsqr, q) + exponent * M_LN2;
z = znum/zden; w = z * z;
x = z + z * w * (POLYNOM2(w,a)/POLYNOM3(w,b));
z = exponent;
x += z * (-2.121944400546905827679e-4);
return x + z * 0.693359375;
}

10
lang/cem/libcc.ansi/math/log10.c
View file

@ -6,16 +6,20 @@
*/
/* $Header$ */
#include <errno.h>
#include <math.h>
#include <errno.h>
#include "localmath.h"
double
log10(double x)
{
if (x <= 0) {
if (x < 0) {
errno = EDOM;
return 0;
return -HUGE_VAL;
}
else if (x == 0) {
errno = ERANGE;
return -HUGE_VAL;
}
return log(x) / M_LN10;

18
lang/cem/libcc.ansi/math/modf.e
View file

@ -5,20 +5,20 @@
*/
/* $Header$ */
mes 2,EM_WSIZE,EM_PSIZE
mes 2,_EM_WSIZE,_EM_PSIZE
#ifndef NOFLOAT
exp $modf
pro $modf,0
lal 0
loi EM_DSIZE
loi _EM_DSIZE
loc 1
loc EM_WSIZE
loc EM_DSIZE
loc _EM_WSIZE
loc _EM_DSIZE
cif
fif EM_DSIZE
lal EM_DSIZE
loi EM_PSIZE
sti EM_DSIZE
ret EM_DSIZE
fif _EM_DSIZE
lal _EM_DSIZE
loi _EM_PSIZE
sti _EM_DSIZE
ret _EM_DSIZE
end
#endif

31
lang/cem/libcc.ansi/math/pow.c
View file

@ -6,13 +6,21 @@
*/
/* $Header$ */
#include <errno.h>
#include <math.h>
#include <float.h>
#include <errno.h>
#include "localmath.h"
double
pow(double x, double y)
{
/* Simple version for now. The Cody and Waite book has
a very complicated, much more precise version, but
this version has machine-dependent arrays A1 and A2,
and I don't know yet how to solve this ???
*/
double dummy;
int result_neg = 0;
if ((x == 0 && y == 0) ||
(x < 0 && modf(y, &dummy) != 0)) {
@ -23,13 +31,26 @@ pow(double x, double y)
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;
result_neg = 1;
}
return val;
x = -x;
}
x = log(x);
if (x < 0) {
x = -x;
y = -y;
}
if (y > M_LN_MAX_D/x) {
errno = ERANGE;
return 0;
}
if (y < M_LN_MIN_D/x) {
errno = ERANGE;
return 0;
}
return exp(log(x) * y);
x = exp(x * y);
return result_neg ? -x : x;
}

119
lang/cem/libcc.ansi/math/sin.c
View file

@ -6,93 +6,76 @@
*/
/* $Header$ */
#include <errno.h>
#include <math.h>
#include "localmath.h"
static double
sinus(double x, int quadrant)
sinus(double x, int cos_flag)
{
/* sin(0.5*pi*x) = x * P(x*x)/Q(x*x) for x in [0,1] */
/* Hart & Cheney # 3374 */
/* Algorithm and coefficients from:
"Software manual for the elementary functions"
by W.J. Cody and W. Waite, Prentice-Hall, 1980
*/
static double p[6] = {
0.4857791909822798473837058825e+10,
-0.1808816670894030772075877725e+10,
0.1724314784722489597789244188e+09,
-0.6351331748520454245913645971e+07,
0.1002087631419532326179108883e+06,
-0.5830988897678192576148973679e+03
static double r[] = {
-0.16666666666666665052e+0,
0.83333333333331650314e-2,
-0.19841269841201840457e-3,
0.27557319210152756119e-5,
-0.25052106798274584544e-7,
0.16058936490371589114e-9,
-0.76429178068910467734e-12,
0.27204790957888846175e-14
};
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;
double xsqr;
double y;
int neg = 0;
if (x < 0) {
quadrant += 2;
x = -x;
neg = 1;
}
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 (cos_flag) {
neg = 0;
y = M_PI_2 + x;
}
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;
else y = x;
/* ??? avoid loss of significance, if y is too large, error ??? */
y = y * M_1_PI + 0.5;
/* Use extended precision to calculate reduced argument.
Here we used 12 bits of the mantissa for a1.
Also split x in integer part x1 and fraction part x2.
*/
#define A1 3.1416015625
#define A2 -8.908910206761537356617e-6
{
double x1, x2;
modf(y, &y);
if (modf(0.5*y, &x1)) neg = !neg;
if (cos_flag) y -= 0.5;
x2 = modf(x, &x1);
x = x1 - y * A1;
x += x2;
x -= n * A2;
x -= y * 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) {
if (x < 0) {
neg = !neg;
x = -x;
}
xsqr = x * x;
x = x * POLYNOM5(xsqr, p) / POLYNOM5(xsqr, q);
return x;
/* ??? avoid underflow ??? */
y = x * x;
x += x * y * POLYNOM7(y, r);
return neg ? -x : x;
}
double

85
lang/cem/libcc.ansi/math/sinh.c
View file

@ -6,33 +6,72 @@
*/
/* $Header$ */
#include <errno.h>
#include <math.h>
#include <float.h>
#include <errno.h>
#include "localmath.h"
static double
sinh_cosh(double x, int cosh_flag)
{
/* Algorithm and coefficients from:
"Software manual for the elementary functions"
by W.J. Cody and W. Waite, Prentice-Hall, 1980
*/
static double p[] = {
-0.35181283430177117881e+6,
-0.11563521196851768270e+5,
-0.16375798202630751372e+3,
-0.78966127417357099479e+0
};
static double q[] = {
-0.21108770058106271242e+7,
0.36162723109421836460e+5,
-0.27773523119650701167e+3,
1.0
};
int negative = x < 0;
double y = negative ? -x : x;
if (! cosh_flag && y <= 1.0) {
/* ??? check for underflow ??? */
y = y * y;
return x + x * y * POLYNOM3(y, p)/POLYNOM3(y,q);
}
if (y >= M_LN_MAX_D) {
/* exp(y) would cause overflow */
#define LNV 0.69316101074218750000e+0
#define VD2M1 0.52820835025874852469e-4
double w = y - LNV;
if (w < M_LN_MAX_D+M_LN2-LNV) {
x = exp(w);
x += VD2M1 * x;
}
else {
errno = ERANGE;
x = HUGE_VAL;
}
}
else {
double z = exp(y);
x = 0.5 * (z + (cosh_flag ? 1.0 : -1.0)/z);
}
return negative ? -x : x;
}
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;
return sinh_cosh(x, 0);
}
double
cosh(double x)
{
if (x < 0) x = -x;
return sinh_cosh(x, 1);
}

2
lang/cem/libcc.ansi/math/sqrt.c
View file

@ -6,8 +6,8 @@
*/
/* $Header$ */
#include <errno.h>
#include <math.h>
#include <errno.h>
#define NITER 5

140
lang/cem/libcc.ansi/math/tan.c
View file

@ -6,117 +6,63 @@
*/
/* $Header$ */
#include <errno.h>
#include <math.h>
#include "localmath.h"
double
tan(double x)
{
/* First reduce range to [0, pi/4].
Then use approximation tan(x*pi/4) = x * P(x*x)/Q(x*x).
Hart & Cheney # 4288
Use: tan(x) = 1/tan(pi/2 - x)
tan(-x) = -tan(x)
tan(x+k*pi) = tan(x)
/* Algorithm and coefficients from:
"Software manual for the elementary functions"
by W.J. Cody and W. Waite, Prentice-Hall, 1980
*/
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;
double y;
static double p[] = {
1.0,
-0.13338350006421960681e+0,
0.34248878235890589960e-2,
-0.17861707342254426711e-4
};
static double q[] = {
1.0,
-0.46671683339755294240e+0,
0.25663832289440112864e-1,
-0.31181531907010027307e-3,
0.49819433993786512270e-6
};
if (negative) x = -x;
/* ??? avoid loss of significance, error if x is too large ??? */
/* 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;
y = x * M_2_PI + 0.5;
/* Use extended precision to calculate reduced argument.
Here we used 12 bits of the mantissa for a1.
Also split x in integer part x1 and fraction part x2.
*/
#define A1 1.57080078125
#define A2 -4.454455103380768678308e-6
{
double x1, x2;
modf(y, &y);
if (modf(0.5*y, &x1)) invert = 1;
x2 = modf(x, &x1);
x = x1 - y * 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;
x -= y * A2;
#undef A1
#undef A2
}
return negative ? -tmp : tmp;
/* ??? avoid underflow ??? */
y = x * x;
x += x * y * POLYNOM2(y, p+1);
y = POLYNOM4(y, q);
if (negative) x = -x;
return invert ? -y/x : x/y;
}

40
lang/cem/libcc.ansi/math/tanh.c
View file

@ -6,19 +6,45 @@
*/
/* $Header$ */
#include <errno.h>
#include <float.h>
#include <math.h>
#include "localmath.h"
double
tanh(double x)
{
if (x <= 0.5*M_LN_MIN_D) {
return -1;
}
/* Algorithm and coefficients from:
"Software manual for the elementary functions"
by W.J. Cody and W. Waite, Prentice-Hall, 1980
*/
static double p[] = {
-0.16134119023996228053e+4,
-0.99225929672236083313e+2,
-0.96437492777225469787e+0
};
static double q[] = {
0.48402357071988688686e+4,
0.22337720718962312926e+4,
0.11274474380534949335e+3,
1.0
};
int negative = x < 0;
if (negative) x = -x;
if (x >= 0.5*M_LN_MAX_D) {
return 1;
x = 1.0;
}
x = exp(x + x);
return (x - 1.0)/(x + 1.0);
#define LN3D2 0.54930614433405484570e+0 /* ln(3)/2 */
else if (x > LN3D2) {
x = 0.5 - 1.0/(exp(x+x)+1.0);
x += x;
}
else {
/* ??? avoid underflow ??? */
double g = x*x;
x += x * g * POLYNOM2(g, p)/POLYNOM3(g, q);
}
return negative ? -x : x;
}