Listing 1: Quad-precision math functions
/* fp.h Floating Point C Extensions header */
#ifndef double_t
#define double_t double
#endif
#ifndef float_t
#ifdef __STDC__
#define float_t float
#else
#define float_t double_t
#endif
#endif
#ifndef _WIDEST_NEED_EVAL
#define _WIDEST_NEED_EVAL 0
#endif
#ifndef _MIN_EVAL_FORMAT /* 80x87 and 6888x: 2; PowerPC, PDP: 1 */
#define _MIN_EVAL_FORMAT 0
#endif
/* math.h standard header */
#ifndef _MATH
#define _MATH
/* macros */
/* some constants from traditional <math.h> and XOPEN */
# define M_E 2.7182818284590452354
# define M_LOG2E 1.4426950408889634074
# define M_LOG10E 0.43429448190325182765
# define M_LN2 0.69314718055994530942
# define M_LN2L 0.6931471805599453094172321214581765680754L
# define M_LN10 2.30258509299404568402
# define M_PI 3.14159265358979323846
# define M_PI_2 1.57079632679489661923
# define M_PI_4 0.78539816339744830962
# define M_1_PI 0.31830988618379067154
# define M_2_PI 0.63661977236758134308
# define M_2_SQRTPI 1.12837916709551257390
# define M_SQRT2 1.41421356237309504880
# define M_SQRT1_2 0.70710678118654752440
/* type definitions */
typedef const union {
unsigned long _W[2];
double _D;
} _Dconst;
/* declarations */
double acos(double);
double asin(double);
double atan(double);
double atan2(double, double);
double ceil(double);
double cos(double);
double cosh(double);
double exp(double);
double fabs(double);
double floor(double);
double fmod(double, double);
double frexp(double, int *);
double ldexp(double, int);
double log(double);
double log10(double);
double modf(double, double *);
double pow(double, double);
double sin(double);
double sinh(double);
double sqrt(double);
double tan(double);
double tanh(double);
long double cosl(long double);
long double sinl(long double);
long double sinhl(long double);
long double coshl(long double);
long double tanhl(long double);
long double tanl(long double);
long double acosl(long double);
long double asinl(long double);
long double atan2l(long double,long double);
long double atanl(long double);
long double logl(long double);
long double log10l(long double);
long double powl(long double,long double);
long double sqrtl(long double);
long double expl(long double);
long double _Logl(long double, int);
double _Log(double, int);
double _Sin(double, unsigned int);
extern _Dconst _Hugeval;
/* macro overrides */
#define atan(x) atan2(x,1.)
#ifndef HPQ
typedef struct{long double s;long double c;long double t;}_scosl;
_scosl _Tanl(long double);
#define tanl(x) _Tanl(x).t
#ifndef SCT
long double _Sinl(long double, unsigned int);
#define sinl(x) _Sinl(x, 0)
#define cosl(x) _Sinl(x, 1)
#else /* for use when 2 or more functions needed */
#define cosl(x) _Tanl(x).c
#define sinl(x) _Tanl(x).s
#endif
#define cos(x) _Sin(x, 1)
#define sin(x) _Sin(x, 0)
#define logl(x) _Logl(x,0)
#define log10(x) _Log(x,1)
#define log(x) _Log(x,0)
#define log10l(x) _Logl(x,1)
#else
long double FTN_QATAN(long double);
long double FTN_QASIN(long double);
long double FTN_QACOS(long double);
long double FTN_QATAN2(long double,long double);
#define atanl(x) FTN_QATAN(x)
#define asinl(x) FTN_QASIN(x)
#define acosl(x) FTN_QACOS(x)
#define atan2l(x,y) FTN_QATAN2(x,y)
long double FTN_QEXP(long double);
long double FTN_QLOG(long double);
long double FTN_QLOG10(long double);
#define expl(x) FTN_QEXP(x)
#define logl(x) FTN_QLOG(x)
#define log10l(x) FTN_QLOG10(x)
long double FTN_QSINH(long double);
long double FTN_QTANH(long double);
long double FTN_QCOSH(long double);
#define sinhl(x) FTN_QSINH(x)
#define coshl(x) FTN_QCOSH(x)
#define tanhl(x) FTN_QTANH(x)
long double FTN_QSIN(long double);
long double FTN_QCOS(long double);
long double FTN_QTAN(long double);
#define sinl(x) FTN_QSIN(x)
#define cosl(x) FTN_QCOS(x)
#define tanl(x) FTN_QTAN(x)
long double FTN_QSQRT(long double);
#define sqrtl(x) FTN_QSQRT(x)
#endif
float acosf(float);
float asinf(float);
float atanf(float);
float atan2f(float, float);
float ceilf(float);
float cosf(float);
float coshf(float);
float expf(float);
float fabsf(float);
float floorf(float);
float fmodf(float, float);
float frexpf(float, int *);
float ldexpf(float, int);
float logf(float);
float log10f(float);
float modff(float, float *);
float powf(float, float);
float sinf(float);
float sinhf(float);
float sqrtf(float);
float tanf(float);
float tanhf(float);
#define frexpf(x,p) (float)frexp(x,p)
#define ldexpf(x,y) (float)ldexp(x,y)
typedef struct{float s;float c;float t;}_scos;
_scos _Sinf(double);
double _Logf(double);
#ifndef HPQ
#define atanf(x) atan2f(x,1.f)
#define logf(x) _Logf(x)
#define log10f(x) (M_LOG10E*_Logf(x))
#define cosf(x) _Sinf(x).c
#define sinf(x) _Sinf(x).s
#ifndef SCT /* separate tanf() more accurate? */
#define tanf(x) _Sinf(x).t
#endif
#endif
/* following macros are unsafe when x,y have side effects */
#if 0 /* these also affect numerical results and do not
improve timing with the HP compiler */
#define acosf(x) atan2f((float)sqrt((1-(double)(x)*(x))),x)
#define asinf(x) atan2f(x,(float)sqrt((1-(double)(x)*(x))))
#define asin(x) atan2(x,sqrt((1-(x))*(1+(x))))
#define acos(x) atan2(sqrt((1-(x))*(1+(x))),x)
#define asinl(x) atan2l(x,sqrtl((1-(x))*(1+(x))))
#define acosl(x) atan2l(sqrtl((1-(x))*(1+(x))),x)
#endif
#define ceilf(x) (float)(long
int)((x)>0?1-.25*FLT_EPSILON+(x):(x))
#if 1
#define fabsf(x) ((x)<0?-(x):(x)) /* not if compiler has own
code */
#define fabs(x) ((x)<0?-(x):(x)) /* not if compiler has own
code */
#endif
#define floorf(x) (float)(long
int)((x)<0?.25*FLT_EPSILON-1+(x):(x))
#define fmodf(x,y) (float)((x)-(long int)((x)/(y))*(double)(y))
#define modff(x,p) ((x)-(*p=(long int)(x)))
#endif
/* xmath.h header */
#include "math.h"
#include <errno.h>
#include <float.h>
#include "fp.h"
typedef union { /* VAX compatibility discarded for simplicity */
double _D;
unsigned long _W[2];
} _Dvar;
/* this section belongs in yvals.h according to Plauger */
#define _D0 0 /* big-endian order */
#define _DOFF 20 /* IEEE double format, using union of longs */
#define _DSIGN 0x80000000 /* not given in article: sign bit
mask */
#define _DBIAS 1022 /* not given in article: exponent field of
0.0 */
#define _FOFF (FLT_MANT_DIG - 1)
/* end of section from yvals.h */
#define _DFRAC ((1<<_DOFF)-1)
#define _DMASK ( ~_DFRAC) /* wrong in article */
#define _DMAX ((1<<(31-_DOFF))-1) /* 15 for short */
#define _DNAN (0x80000000|_DMAX<<_DOFF|1<<(_DOFF-1)) /* 8000 */
#define _D1 (1-_D0) /* word order VAX */
#define FINITE -1
#define INF 1
#define NAN 2
#define NBITS (32+_DOFF)
#define DSIGN(x) (((unsigned long *)&(x))[_D0]&_DSIGN)
#define HUGE_EXP (int)(_DMAX * 900L / 100) /* correct bug in
Plauger */
#define HUGE_RAD 3.14e30
#define SAFE_EXP (_DMAX>>1)
/* declarations */
int _Dint();
int _Dnorm();
int _Dscale();
double _Dtento();
int _Dtest();
#define _Dtest(x)
((x)>DBL_MAX||(x)<-DBL_MAX?INF:(x)!=(x)?NAN:(x)!=0?FINITE:0)
#define _Ltest(x)
((x)>LDBL_MAX||(x)<-LDBL_MAX?INF:(x)!=(x)?NAN:(x)!=0?FINITE:0)
int _Dunscale();
int _Exp();
int _Ldunscale();
extern _Dconst _Hugeval;
/* atanl() */
/* Copyright (c) 1985 Regents of the University of California.
Use and reproduction of this software are granted in accordance
with
the terms and conditions specified in the Berkeley Software
License
Agreement (in particular, this entails acknowledgement of the
programs' source, and inclusion of this notice) with the
additional
understanding that all recipients should regard themselves as
participants in an ongoing research project and hence should
feel
obligated to report their experiences (good or bad) with these
elementary function codes, using "sendbug 4bsd-bugs@BERKELEY",
to the
authors.
*/
#ifndef lint
static char sccsid[] = "@(#)atan2l.c (T Prince) 6/3/95";
#endif /* not lint */
/* ATAN2(Y,X) RETURN ARG (X+iY)
* DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
* CODED IN C BY K.C. NG, 1/8/85;
* REVISED BY K.C. NG on 2/7/85, 2/13/85, 3/7/85, 3/30/85,
6/29/85.
* long double by T C Prince 6/95
*
* Method :
1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
2. Reduce x to positive by (if x and y are unexceptional):
ARG (x+iy) = arctan(y/x) ... if x > 0,
ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
3. the argument
is further reduced to one of the following intervals and
the
arctangent of y/x is evaluated by the corresponding
formula:
[0,1/3] atan(y/x) = t -
t^3*(a1+t^2*(a2+...a10+t^2*a11)
[1/3,11/16] atan(y/x) = atan(1/2) + atan(
(y-x/2)/(x+y/2) )
[11/16,19/16] atan(y/x) = atan( 1 ) + atan( (y-x)/(x+y)
)
[19/16,3] atan(y/x) = atan(3/2) +
atan((y-1.5x)/(x+1.5y))
[3,INF] atan(y/x) = atan(INF) + atan( -x/y )
The decimal values may be used, provided that the compiler will
convert from decimal to binary accurately enough.
*/
#define ath .463647609000806116214256231461214402028537L
#define at1 .982793723247329067985710611014666014496877L
#define athfhi (double)ath
#include <float.h>
#define FUDGE (1+LDBL_EPSILON/DBL_EPSILON)
#define athflo (ath - athfhi)*FUDGE
#define at1fhi (double)at1
#define at1flo (at1-at1fhi)*FUDGE
#define PI 3.141592653589793238462643383279502884197169L
#include "math.h"
#include <assert.h>
#include "fp.h"
long double atan2l(long double y, long double x)
{
long double t, z, signy, hi, lo;
double_t a11;
int signx;
/* Copy down the sign of y and x */
signy = 1.;
if (y < 0) {
y = -y;
signy = -1.;
}
if (signx = x < 0)
x = -x;
if ((a11 = t = y / x) < 3) {
/* In case T is in [0,1/3] */
hi = 0;
lo = 0;
if (a11 > 1 / 3.) {
if (a11 < 11 / 16.) { /* T is in [1/3,11/16] */
t = (y + y - x) / (x + x + y);
hi = athfhi;
lo = athflo;
} else {
long double one5 = 1.5;
if (a11 < 19 / 16.) {
/* T is in [11/16,19/16] */
t = (y - x) / (x + y);
hi = (double) PI *.25;
lo = (PI * .25 - (double) PI * .25) * FUDGE;
} else { /* T is in [19/16,3] */
t = (y - x * one5) / (x + y * one5);
hi = at1fhi;
lo = at1flo;
}
}
}
} else { /* T >= 3 */
/* 0/0 produces NaN (agrees with IEEE, not UNIX tradition) */
t = -x / y;
hi = (double) PI *.5;
lo = (PI * .5 - (double) PI * .5) * FUDGE;
}
#ifdef DEBUG
assert(t <= 1. / 3 || !(t == t));
#endif
a11 = z = t * t;
z = (t + (lo - t * z *
(.333333333333333333333333333333326238L +
z * (-.199999999999999999999999999989692064L +
z * (.142857142857142857142857136902464695L
+
z * (-.111111111111111111111109286933220956L
+
z * (.090909090909090909090566375305311242L
+
z * (-.076923076923076923033848105066780741L
+
z * (.066666666666666662832160589417621446L
+
z * (-.058823529411764454382902973039123863L
+
z * (.052631578947355920010163119267195941L
+
z * (-.047619047618567045788940099947698386L
+
z * (.043478260855068916374201546117917056L
+
z * (-.039999999653506693442185507209875834L
+
z * (.037037030436227061574079324145862133L
+
z * (-.034482658185047909429270585023057681L
+
z * (.032256845866214999368038205000384542L
+
z * (-.030291305255230381718828679590094691L
+
z * (.028482910936322424263248147209081714L
+
z * (-.026511385701218695500640341310485698L
+
a11 * (.023381127183313153463575187882151062
+
a11 *
(-.017211039939973303258669077432249085 +
a11 * .007572335682370299569957516923702074
)))))))))))))))))))))) + hi;
return (signx ? PI - z : z) * signy;
}
long double (atanl) (long double x) {
return atan2l(x, 1.);
}
long double (asinl) (long double x) {
return atan2l(x, sqrtl((1 - x) * (1 + x)));
}
long double (acosl) (long double x) {
return atan2l(sqrtl((1 - x) * (1 + x)), x);
}
/* coshl() */
#include "math.h"
#include <float.h>
#include "fp.h"
long double (coshl) (long double x) {
double_t xd = x;
if (xd < 0)
x = -x;
if (xd > 11355) /* avoid premature overflow */
return expl(x - .6931471805599453094172321214581765680754L);
x = expl(x); /* sometimes more accurate */
return .5 * x + 0.5 / x;
}
/* expl() */
long double expl(long double x)
{
#define max(i,j) ((i)>(j)?(i):(j))
#define min(i,j) ((i)<(j)?(i):(j))
#include <float.h>
#include <errno.h>
#include <limits.h>
#include <assert.h>
#include "fp.h"
#define ln2 .6931471805599453094172321214581765680754L
#define ln2f (float)ln2
#define LOFF 16
int mi, msign;
double_t g = x;
long double x2;
if (g != g) {
errno = EDOM;
return x;
}
if (g > INT_MAX / 2)
return LDBL_MAX * 2; /* (long double)HUGE_VAL */
if (g < -INT_MAX / 2)
return 0;
g = mi = 1 / ln2 * x + (g < 0 ? -.5 : .5);
/* ln2f*g is exact without long double */
msign = (g = x = (x - ln2f * g) - g * (ln2 - ln2f)) < 0;
#ifdef DEBUG
assert(g >= -.5 * ln2f && g <= .5 * ln2f);
#endif
#ifdef RAT
/* approximation of less than 21 digits precision; this form
of approximation should be more efficient than Chebyshev
economized
polynomial */
x2 = x * x;
x += x * x2 * (.030302142733924400603L + x2 *
.000126241493006006923L);
/* normal code reads .5 + x / ((x2 ..- x) + 2); with mi and
the limits
* incremented */
x = 1 + 1 / ((x2 * (.227270952134515470443L + x2 *
(.002525062330554872359L + x2 * .000003005196267820965)) -
x) + 2)
* (x + x);
#else
x = 1 + (x + x * x * (.5 + x *
(.1666666666666666666666666666666666L +
x * (.0416666666666666666666666666666879L +
x * (.0083333333333333333333333333333364L +
x * (.0013888888888888888888888888861227L +
x * (.0001984126984126984126984126981051L +
x * (.0000248015873015873015873017716641L +
x * (.0000027557319223985890652557486416L +
x * (.0000002755731922398589065184059538L +
x * (.0000000250521083854417187744983345L +
x * (.0000000020876756987868100716035797L +
x * (.0000000001605904383682161577574625L +
x * (.0000000000114707455977270220438353L +
x * (.0000000000007647163731818179510622L +
x * (.0000000000000477947733508877057479L +
g * (.0000000000000028114572558302611040 +
g * .0000000000000001561919020418582642 +
g * g * (.0000000000000000082206267924090703
+
g * (.0000000000000000004116196036662569 +
g *
.0000000000000000000196003744570299)))))))))))))))))));
#endif
*(int *) &x += (max(LDBL_MIN_EXP, min(LDBL_MAX_EXP + msign -
1,
mi)) << LOFF);
/* x.ex+=mi; with limiting to prevent exponent wraparound */
if (mi < LDBL_MIN_EXP) {
/* complete underflow by multiplication */
long double g = 1.;
*(int *) &g += (max(LDBL_MIN_EXP,
mi - LDBL_MIN_EXP)) << LOFF;
x *= g;
}
if (mi > LDBL_MAX_EXP + msign - 1) {
x = LDBL_MAX * 2;
errno = ERANGE;
}
return x;
}
/* logl function */
#include <float.h>
#include "xmath.h"
#include <errno.h>
long double (_Logl) (long double x, int dec) { /* compute
ln(x) */
int xexp, inc;
long double z, w, x1;
double_t x4;
const static long double y = 0.70710678118654752440;
long x0, y0 = ((long *) &y)[_D0];
#define LBIAS (LDBL_MAX_EXP-2)
#define LOFF 16
#define LMASK ((-1)<<LOFF)
#define log2
.693147180559945309417232121458176568075500134360255254120679L
#define lg10 .4342944819032518276511289189166050822943970058036L
#define l2d (double)log2
#define l210d (double)(log2*lg10)
#include <assert.h>
if (x >= 0) { /* ln(negative) undefined */
xexp = 0;
if (x <= LDBL_MAX) {
if (x < LDBL_MIN) { /* take care of subnormal */
x *= 2 / LDBL_EPSILON;
xexp = -LDBL_MANT_DIG;
}
/* extract exponent field */
xexp += ((((unsigned int *) &x)[_D0] & LMASK) >> LOFF) -
LBIAS;
/* adjust by 0 or 1 according to mantissa range */
xexp -= inc = (unsigned long) ((x0 = ((unsigned int *)
&x)[_D0]
& ~LMASK | (LBIAS << LOFF)) - y0) >> 31;
/* set in reduced range */
((unsigned int *) &x)[_D0] = x0 + (inc << LOFF);
/* If x < sqrt(1/2) shift up by factor 2 */
#ifdef DEBUG
assert(x >= y && x <= y * 2);
#endif
/* Represent as log(2^xexp*x), sqrt(.5)<= x <= sqrt(2)
*/
z = (x1 = x - 1) / (x + 1);
x4 = w = z * z;
/* represent 2z as x1-z*x1 to reduce error */
z *= x1 - w * (.6666666666666666666666666666666651039L +
w * (.4000000000000000000000000000050915483L +
w * (.2857142857142857142857142791488061900L
+
w * (.2222222222222222222222266773772078210L
+
w * (.1818181818181818181799844854335815866L
+
w * (.1538461538461538466513537079051311735L
+
w * (.1333333333333332395423334986652396869L
+
w * (.1176470588235421264298007870862879176L
+
w * (.1052631578934701082005085235538617685L
+
w * (.0952380953319759987556163464587572459L
+
w * (.0869565165443288621007967117121022632L
+
w * (.0800002135279663799757893380567682226L
+
x4 * (.0740676564232376969143453696329673648
+
x4 * .0691021859333692522131244214299015772
+
x4 * x4 * (.0625725802128490075702915695648370866
+
x4 *
.0769349179526519185536102742270211149))))))))))))));
w = xexp;
return dec ? (lg10 * x1 - lg10 * z + (log2 * lg10 -
l210d) * w)
+ l210d * w : (x1 - z + w * (log2 - l2d)) + l2d * w;
}
errno = ERANGE;
return x; /* Inf */
}
errno = EDOM; /* NaN */
return x != x ? x : (LDBL_MAX * 2) * 0; /* NaN */
}
long double (logl) (long double x) {
return _Logl(x, 0);
}
long double (log10l) (long double x) {
return _Logl(x, 1);
}
/* powl function */
#include "xmath.h"
#include "math.h"
#include "fp.h"
#include <limits.h>
#include <float.h>
#define LBIAS (LDBL_MAX_EXP-2)
#define LOFF 16
#define LMASK ((-1)<<LOFF)
#define log2 .693147180559945309417232121458176568075500134L
long double (powl) (long double x, long double y) { /* compute
x^y */
double_t yd = y;
if (yd <= INT_MAX && yd >= INT_MIN) {
int exp = yd, xexp, neg, rexp;
long double z, yi = x, yf = y - exp, w;
/* Multiply out to allow x<=0 or y==0 */
register int expr = exp;
if (x == 0 && yd > 0)
return x;
if (neg = (expr < 0)) {
expr = -expr;
}
z = 1.; /* powl(0,0) = 1 */
if ((rexp = (yf == 0)) && neg)
yi = 1 / yi;
if (expr & 1)
z = yi;
for (; expr > 1; z *= yi)
do
yi *= yi;
while (((expr >>= 1) & 1) == 0);
if (rexp)
return z;
/* -Inf falls through to NaN, what about +Inf ? */
if (x > 0) {
xexp = 0;
if (x >= LDBL_MIN) {
/* no danger of underflow, so
* calculate log(x) as log(2^n*xn) */
xexp = (rexp = ((((unsigned int *) &x)[_D0] &
LMASK) >> LOFF) - LBIAS);
((unsigned int *) &x)[_D0] -= rexp << LOFF;
}
w = xexp * yf;
x = expl(logl(x) * yf + (w - (xexp = w)) * log2);
((unsigned int *) &x)[_D0] += xexp << LOFF;
return neg ? x / z : z * x;
}
errno = EDOM;
return x != x ? x : (LDBL_MAX * 2) * 0;
}
return expl(logl(x) * y);
}
/* sinhl function */
#include "math.h"
#include "xmath.h"
#include "fp.h"
long double (sinhl) (long double x) { /* compute sinhl(x) */
int neg;
double_t xd = x;
switch (_Dtest(xd)) { /* test for special codes */
case NAN:
errno = EDOM;
return (x);
case INF:
errno = ERANGE;
return (x);
default: /* finite */
{ /* compute sinh(finite) */
if (neg = (xd < 0.0)) {
x = -x;
xd = -xd;
}
if (xd < 1) { /* |x| < 1 */
long double y = x * x;
double_t yd = y;
x += x * y * (.16666666666666666666666666666666573L +
y * (.00833333333333333333333333333338555L +
y * (.00019841269841269841269841269726379L +
y * (.00000275573192239858906525574505191L +
y * (.00000002505210838544171877496283537L +
y * (.00000000016059043836821614638343470L +
y * (.00000000000076471637318198050919003L +
y * (.00000000000000281145725434779709982L +
y * (.00000000000000000822063524350084984L +
yd * (.00000000000000000001957294395545097 +
yd * (.00000000000000000000003867997525529 +
yd *
.00000000000000000000000006506911776)))))))))));
} else if (xd > 11355) /* avoid premature overflow */
x = expl(x -
.6931471805599453094172321214581765680754L);
else {
x = expl(x); /* sometimes more accurate */
x = .5 * x - 0.5 / x;
}
return (neg ? -x : x);
}
}
}
/* tanhf function */
#include "xmath.h"
float (tanhf) (float x) {
/* Fetch some constants early to get pipelines started */
double_t xx;
float_t minthreshold = -1;
if (x <= 1) {
/* Postpone 2nd conditional so it can run parallel */
#define c4 15.22528717
#define c6 430.2689631
float_t c2 = c6 / c4;
double_t x2 = x * x, c3 = (107.2745303 - c6) /
c4, c5 = (1 - c4) / c4;
if (x >= minthreshold) return x + x * x2 * (c3 + x2 * c5) /
(968.9833252 / c4 + x2 * (c2 + x2));
}
xx = x + x;
if (x == x) {
if (_Expf(&xx) < DBL_MANT_DIG) /* Promote final calc to
double */
return(xx - 1) / (xx + 1);
return 1; /* large x */
}
errno = EDOM; /* NaN */
return x;
}
/* tanhl function */
#include "math.h"
#include "xmath.h"
#include "fp.h"
/* coefficients, after Cody & Waite, Chapter 13 */
long double (tanhl) (long double x) { /* compute tanh(x) */
double_t gd = x, ln3by2 = 0.54930614433405484570;
switch (_Dtest(gd)) { /* test for special codes */
case NAN:
errno = EDOM;
return x;
case INF:
return (gd < 0 ? -1 : 1);
default: /* compute tanhl(finite) */
if (gd < ln3by2 && gd > -ln3by2) { /* |x| < ln(3)/2 */
/* rational approximation for non-linear term should be
more
* economical */
long double g = x * x;
gd = g;
#ifdef RAT
/* this rational approx good only for 19 digits */
#define th06 28.07024027810745285L
#define th0 135867.2643694332105L/(1-th06)
#define th1 (135867.2643694332105L/(1-th06)-th0)
#define th2 62702.82599031282509L/(1-th06)
#define th3 (17413.73786716842593L/(1-th06)-th2)
#define th4 3164.784028589226048L/(1-th06)
#define th5 (379.4772810758125617L/(1-th06)-th4)
#define th6 (double)(28.07024027810745285L/(1-th06))
#define th7 (1.-th6)
return x + x * g * (th3 + g * (th5 + g)) /
(th0 + g * (th2 + g * (th4 + th6 * gd)));
#else
return x + x * g *
(-.33333333333333333333333333333331236L +
g * (.13333333333333333333333333332106849L +
g * (-.05396825396825396825396825111388091L
+
g * (.02186948853615520282186913593999983L +
g * (-.00886323552990219656883648897245587L
+
g * (.00359212803657248101556579466968290L +
g * (-.00145583438705131821918764815034089L
+
g * (.00059002744094558467345621705938886L +
g * (-.00023912911424352596903127485564400L
+
g * (.00009691537956887717041539595368586L +
g * (-.00003927832387813678621340928324239L
+
g * (.00001591890501807617514370588356740L +
g * (-.00000645168880852671761352918937542L
+
g * (.00000261476854493288766376082416958L +
g * (-.00000105971337648849907461636231708L
+
g * (.00000042943523810581307720525608556L +
g * (-.00000017388079025122267319790753556L
+
g * (.00000007005993273168794852063765611L +
gd * (-.00000002759864009681181986545670195
+
gd * .00000001004869246620310962383899486 +
gd * gd * (-.00000000294304745646220434203258075
+
gd *
.00000000049586245014711828294455498))))))))))))))))))));
#endif
}
return 1 - 2 / (expl(x + x) + 1);
}
}
/* tanl function */
#include "xmath.h"
#include "fp.h"
#include <float.h>
#include <assert.h>
_scosl _Tanl(long double x)
{ /* compute tan(x) */
long double g, gd, y;
long quad;
int qt;
double_t gr = 1 / LDBL_EPSILON, grd;
_scosl res;
g = x * (2 / 3.141592653589793238462643383279502884L);
if (FLT_ROUNDS > 0) {
if (x < 0)
gr = -gr;
if (FLT_ROUNDS == 1) /* ANINT(x*2/pi) */
g = (g + gr) - gr;
else
g = ((g + (FLT_ROUNDS == 2 ?
-.5 : .5)) + gr) - gr;
} else { /* Assume FLT_ROUNDS ==0 for positive
numbers */
g = x < 0 ?
gr - ((.5 - g) + gr) : ((g + .5) + gr) - gr;
}
gd = gr * 4;
quad = g - ((g + gd) - gd); /* prevent overflow */
#ifdef DEBUG
assert(abs(quad) <= 2);
#endif
#define Pi2 (.5* 3.141592653589793238462643383279502884L)
g = (x - g * (double) Pi2) - g * (Pi2 - (double) Pi2);
gr = y = g * g;
/* copied from sinl(), needs further economization if not
using sinl(),
* cosl() results */
#define c2 -.49999999999999999999999999999999992447L
#define c4 .04166666666666666666666666666665981199L
#define c6 -.00138888888888888888888888888864445448L
#define c8 .00002480158730158730158730158277331056L
#define c10 -.00000027557319223985890652552327694548L
#define c12 .00000000208767569878680989756792033298L
#define c14 -.00000000001147074559772972304069435075L
#define c16 .00000000000004779477332386842803248038L
#define c18 -.00000000000000015619206967379035718829L
#define c20 .00000000000000000041103174419983283832L
#define c22 -.00000000000000000000088966157295338194
#define c24 .00000000000000000000000160182379249085
#define c3 -.16666666666666666666666666666665583L
#define c5 .008333333333333333333333333332496349L
#define c7 -.00019841269841269841269841267308313L
#define c9 .000002755731922398589065255335933933L
#define c11 -.00000002505210838544171877139980846L
#define c13 .000000000160590438368216124640564804L
#define c15 -.00000000000076471637318189946861586L
#define c17 .000000000000002811457254134541453643L
#define c19 -.00000000000000000822063488884246675
#define c21 .000000000000000000019572556231225041
#define c23 .00000000000000000000003844380218120
grd = gr * (c19 + gr * (c21 - gr * c23));
gd = y * (c2 + y * (c4 + y * (c6 + y * (c8 + y * (c10 + y *
(c12 +
y * (c14 + y * (c16 + y * (c18 + y * (c20 +
gr * (c22 + gr * c24)))))))))));
y *= (c3 + y * (c5 + y * (c7 + y * (c9 +
y * (c11 + y * (c13 + y * (c15 + y * (c17 +
grd))))))));
if ((unsigned int) quad & 1) { /* reduced arg near +- Pi/2 */
res.s = -1 - gd; /* sin (near -1) */
res.c = g + g * y; /* cos */
res.t = res.s / res.c; /* tan */
if (quad == 1) { /* g near Pi/2 */
res.s = -res.s;
res.c = -res.c;
}
} else { /* reduced arg near +- Pi or 0 */
res.c = 1 + gd; /* cos (near 1) */
res.s = g + g * y; /* sin */
res.t = g + g * (y - gd) / res.c; /* tan */
if (quad & 2) { /* reduced arg near +-Pi */
res.s = -res.s;
res.c = -res.c;
}
}
return res;
}
#if 0 /* don't use unless wanting 2 functions
*/
long double (cosl) (long double x) {
return _Tanl(x).c;
}
long double (sinl) (long double x) {
return _Tanl(x).s;
}
#endif
long double (tanl) (long double x) {
return _Tanl(x).t;
}
/* _Sinl function */
#include "math.h"
#include <float.h>
#include "xmath.h"
long double _Sinl(long double x, unsigned int qoff)
{ /* sin(x) or cos(x) */
#define pi 3.141592653589793238462643383279502884L
long double g = 2 / pi * x, g2;
double_t gr = 1 / LDBL_EPSILON, g2d;
int quad;
#ifndef HPQ
if (FLT_ROUNDS > 0) {
if (x < 0)
gr = -gr;
if (FLT_ROUNDS == 1) /* ANINT(x*2/pi) */
g = (g + gr) - gr;
else
g = ((g + (FLT_ROUNDS == 2 ?
-.5 : .5)) + gr) - gr;
} else { /* Assume FLT_ROUNDS ==0 for positive
numbers */
g = x < 0 ?
gr - ((.5 - g) + gr) : ((g + .5) + gr) - gr;
}
g2 = gr * 4;
quad = g - ((g + g2) - g2); /* prevent overflow */
#else /* HP FORTRAN math library functions */
long double FTN_QNINT(long double),FTN_QMOD(long double,long
double);
quad = FTN_QMOD(g = FTN_QNINT(g),4); /*
g=ANINT(g);quad=MOD(g,4Q0) */
#endif
/* Get remainder after subtracting nearest multiple of pi/2
*/
#define hpi (.5 * pi)
#define hpid (double)hpi
g = (x - hpid * g) - g * (hpi - hpid);
g2d = g2 = g * g;
if ((qoff += quad) & 1) /* cosine series */
g = 1 + g2 * (-.49999999999999999999999999999999992447L +
g2 * (.04166666666666666666666666666665981199L +
g2 *
(-.00138888888888888888888888888864445448L +
g2 * (.00002480158730158730158730158277331056L
+
g2 * (-.00000027557319223985890652552327694548L
+
g2 * (.00000000208767569878680989756792033298L
+
g2 * (-.00000000001147074559772972304069435075L
+
g2 * (.00000000000004779477332386842803248038L
+
g2 * (-.00000000000000015619206967379035718829L
+
g2 * (.00000000000000000041103174419983283832L
+
g2d * (-.00000000000000000000088966157295338194
+
g2d *
.00000000000000000000000160182379249085)))))))))));
else /* sine series */
g += g * g2 * (-.16666666666666666666666666666665583L +
g2 * (.008333333333333333333333333332496349L +
g2 * (-.00019841269841269841269841267308313L
+
g2 * (.000002755731922398589065255335933933L
+
g2 * (-.00000002505210838544171877139980846L
+
g2 * (.000000000160590438368216124640564804L
+
g2 * (-.00000000000076471637318189946861586L
+
g2 * (.000000000000002811457254134541453643L
+
g2d * (-.00000000000000000822063488884246675
+
g2d * (.000000000000000000019572556231225041
-
g2d *
.00000000000000000000003844380218120))))))))));
return qoff & 2 ? -g : g;
}
long double (sinl) (long double x) {
return _Sinl(x, 0);
}
long double (cosl) (long double x) {
return _Sinl(x, 1);
}