jyapayne/musl
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

math: sinh.c cleanup similar to the cosh one

Browse source 
comments are kept in the double version of the function
This commit is contained in:
Szabolcs Nagy 2012-12-16 19:49:55 +01:00
commit f143458223
3 changed files with 72 additions and 171 deletions
Download patch file Download diff file Expand all files Collapse all files

90
src/math/sinh.c
View file

@ -1,71 +1,39 @@
/* origin: FreeBSD /usr/src/lib/msun/src/e_sinh.c */
/*
* ====================================================
* 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.
* ====================================================
*/
/* 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 "libm.h"
static const double huge = 1.0e307;
/* sinh(x) = (exp(x) - 1/exp(x))/2
* = (exp(x)-1 + (exp(x)-1)/exp(x))/2
* = x + x^3/6 + o(x^5)
*/
double sinh(double x)
{
double t, h;
int32_t ix, jx;
/* High word of |x|. */
GET_HIGH_WORD(jx, x);
ix = jx & 0x7fffffff;
/* x is INF or NaN */
if (ix >= 0x7ff00000)
return x + x;
union {double f; uint64_t i;} u = {.f = x};
uint32_t w;
double t, h, absx;
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 */
/* raise inexact, return x */
if (huge+x > 1.0)
if (u.i >> 63)
h = -h;
/* |x| */
u.i &= (uint64_t)-1/2;
absx = u.f;
w = u.i >> 32;
/* |x| < log(DBL_MAX) */
if (w < 0x40862e42) {
t = expm1(absx);
if (w < 0x3ff00000) {
if (w < 0x3ff00000 - (26<<20))
/* note: inexact is raised by expm1 */
/* note: this branch avoids underflow */
return x;
t = expm1(fabs(x));
if (ix < 0x3ff00000)
return h*(2.0*t - t*t/(t+1.0));
return h*(t + t/(t+1.0));
return h*(2*t - t*t/(t+1));
}
/* note: |x|>log(0x1p26)+eps could be just h*exp(x) */
return h*(t + t/(t+1));
}
/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
if (ix < 0x40862E42)
return h*exp(fabs(x));
/* |x| in [log(maxdouble), overflowthresold] */
if (ix <= 0x408633CE)
return h * 2.0 * __expo2(fabs(x)); /* h is for sign only */
/* |x| > overflowthresold, sinh(x) overflow */
return x*huge;
/* |x| > log(DBL_MAX) or nan */
/* note: the result is stored to handle overflow */
t = 2*h*__expo2(absx);
return t;
}