Commit f399fc8d authored by Raymond Toy's avatar Raymond Toy
Browse files

Import hyperbolic functions from fdlibm, as is.

parent d29ef8f7
/* @(#)e_cosh.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* __ieee754_cosh(x)
* Method :
* mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
* 1. Replace x by |x| (cosh(x) = cosh(-x)).
* 2.
* [ exp(x) - 1 ]^2
* 0 <= x <= ln2/2 : cosh(x) := 1 + -------------------
* 2*exp(x)
*
* exp(x) + 1/exp(x)
* ln2/2 <= x <= 22 : cosh(x) := -------------------
* 2
* 22 <= x <= lnovft : cosh(x) := exp(x)/2
* lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2)
* ln2ovft < x : cosh(x) := huge*huge (overflow)
*
* Special cases:
* cosh(x) is |x| if x is +INF, -INF, or NaN.
* only cosh(0)=1 is exact for finite x.
*/
#include "fdlibm.h"
#ifdef __STDC__
static const double one = 1.0, half=0.5, huge = 1.0e300;
#else
static double one = 1.0, half=0.5, huge = 1.0e300;
#endif
#ifdef __STDC__
double __ieee754_cosh(double x)
#else
double __ieee754_cosh(x)
double x;
#endif
{
double t,w;
int ix;
unsigned lx;
/* High word of |x|. */
ix = __HI(x);
ix &= 0x7fffffff;
/* x is INF or NaN */
if(ix>=0x7ff00000) return x*x;
/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
if(ix<0x3fd62e43) {
t = expm1(fabs(x));
w = one+t;
if (ix<0x3c800000) return w; /* cosh(tiny) = 1 */
return one+(t*t)/(w+w);
}
/* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
if (ix < 0x40360000) {
t = __ieee754_exp(fabs(x));
return half*t+half/t;
}
/* |x| in [22, log(maxdouble)] return half*exp(|x|) */
if (ix < 0x40862E42) return half*__ieee754_exp(fabs(x));
/* |x| in [log(maxdouble), overflowthresold] */
lx = *( (((*(unsigned*)&one)>>29)) + (unsigned*)&x);
if (ix<0x408633CE ||
(ix==0x408633ce)&&(lx<=(unsigned)0x8fb9f87d)) {
w = __ieee754_exp(half*fabs(x));
t = half*w;
return t*w;
}
/* |x| > overflowthresold, cosh(x) overflow */
return huge*huge;
}
/* @(#)e_sinh.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* __ieee754_sinh(x)
* Method :
* mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
* 1. Replace x by |x| (sinh(-x) = -sinh(x)).
* 2.
* E + E/(E+1)
* 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x)
* 2
*
* 22 <= x <= lnovft : sinh(x) := exp(x)/2
* lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2)
* ln2ovft < x : sinh(x) := x*shuge (overflow)
*
* Special cases:
* sinh(x) is |x| if x is +INF, -INF, or NaN.
* only sinh(0)=0 is exact for finite x.
*/
#include "fdlibm.h"
#ifdef __STDC__
static const double one = 1.0, shuge = 1.0e307;
#else
static double one = 1.0, shuge = 1.0e307;
#endif
#ifdef __STDC__
double __ieee754_sinh(double x)
#else
double __ieee754_sinh(x)
double x;
#endif
{
double t,w,h;
int ix,jx;
unsigned lx;
/* High word of |x|. */
jx = __HI(x);
ix = jx&0x7fffffff;
/* x is INF or NaN */
if(ix>=0x7ff00000) return x+x;
h = 0.5;
if (jx<0) h = -h;
/* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
if (ix < 0x40360000) { /* |x|<22 */
if (ix<0x3e300000) /* |x|<2**-28 */
if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
t = expm1(fabs(x));
if(ix<0x3ff00000) return h*(2.0*t-t*t/(t+one));
return h*(t+t/(t+one));
}
/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
if (ix < 0x40862E42) return h*__ieee754_exp(fabs(x));
/* |x| in [log(maxdouble), overflowthresold] */
lx = *( (((*(unsigned*)&one)>>29)) + (unsigned*)&x);
if (ix<0x408633CE || (ix==0x408633ce)&&(lx<=(unsigned)0x8fb9f87d)) {
w = __ieee754_exp(0.5*fabs(x));
t = h*w;
return t*w;
}
/* |x| > overflowthresold, sinh(x) overflow */
return x*shuge;
}
/* @(#)s_tanh.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* Tanh(x)
* Return the Hyperbolic Tangent of x
*
* Method :
* x -x
* e - e
* 0. tanh(x) is defined to be -----------
* x -x
* e + e
* 1. reduce x to non-negative by tanh(-x) = -tanh(x).
* 2. 0 <= x <= 2**-55 : tanh(x) := x*(one+x)
* -t
* 2**-55 < x <= 1 : tanh(x) := -----; t = expm1(-2x)
* t + 2
* 2
* 1 <= x <= 22.0 : tanh(x) := 1- ----- ; t=expm1(2x)
* t + 2
* 22.0 < x <= INF : tanh(x) := 1.
*
* Special cases:
* tanh(NaN) is NaN;
* only tanh(0)=0 is exact for finite argument.
*/
#include "fdlibm.h"
#ifdef __STDC__
static const double one=1.0, two=2.0, tiny = 1.0e-300;
#else
static double one=1.0, two=2.0, tiny = 1.0e-300;
#endif
#ifdef __STDC__
double tanh(double x)
#else
double tanh(x)
double x;
#endif
{
double t,z;
int jx,ix;
/* High word of |x|. */
jx = __HI(x);
ix = jx&0x7fffffff;
/* x is INF or NaN */
if(ix>=0x7ff00000) {
if (jx>=0) return one/x+one; /* tanh(+-inf)=+-1 */
else return one/x-one; /* tanh(NaN) = NaN */
}
/* |x| < 22 */
if (ix < 0x40360000) { /* |x|<22 */
if (ix<0x3c800000) /* |x|<2**-55 */
return x*(one+x); /* tanh(small) = small */
if (ix>=0x3ff00000) { /* |x|>=1 */
t = expm1(two*fabs(x));
z = one - two/(t+two);
} else {
t = expm1(-two*fabs(x));
z= -t/(t+two);
}
/* |x| > 22, return +-1 */
} else {
z = one - tiny; /* raised inexact flag */
}
return (jx>=0)? z: -z;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment