Back to index

glibc  2.9
e_sinhl.c
Go to the documentation of this file.
00001 /* e_asinhl.c -- long double version of e_asinh.c.
00002  * Conversion to long double by Ulrich Drepper,
00003  * Cygnus Support, drepper@cygnus.com.
00004  */
00005 
00006 /*
00007  * ====================================================
00008  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00009  *
00010  * Developed at SunPro, a Sun Microsystems, Inc. business.
00011  * Permission to use, copy, modify, and distribute this
00012  * software is freely granted, provided that this notice
00013  * is preserved.
00014  * ====================================================
00015  */
00016 
00017 #if defined(LIBM_SCCS) && !defined(lint)
00018 static char rcsid[] = "$NetBSD: $";
00019 #endif
00020 
00021 /* __ieee754_sinhl(x)
00022  * Method :
00023  * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
00024  *     1. Replace x by |x| (sinhl(-x) = -sinhl(x)).
00025  *     2.
00026  *                                                 E + E/(E+1)
00027  *         0        <= x <= 25     :  sinhl(x) := --------------, E=expm1l(x)
00028  *                                                        2
00029  *
00030  *         25       <= x <= lnovft :  sinhl(x) := expl(x)/2
00031  *         lnovft   <= x <= ln2ovft:  sinhl(x) := expl(x/2)/2 * expl(x/2)
00032  *         ln2ovft  <  x        :  sinhl(x) := x*shuge (overflow)
00033  *
00034  * Special cases:
00035  *     sinhl(x) is |x| if x is +INF, -INF, or NaN.
00036  *     only sinhl(0)=0 is exact for finite x.
00037  */
00038 
00039 #include "math.h"
00040 #include "math_private.h"
00041 
00042 #ifdef __STDC__
00043 static const long double one = 1.0, shuge = 1.0e4931L;
00044 #else
00045 static long double one = 1.0, shuge = 1.0e4931L;
00046 #endif
00047 
00048 #ifdef __STDC__
00049        long double __ieee754_sinhl(long double x)
00050 #else
00051        long double __ieee754_sinhl(x)
00052        long double x;
00053 #endif
00054 {
00055        long double t,w,h;
00056        u_int32_t jx,ix,i0,i1;
00057 
00058     /* Words of |x|. */
00059        GET_LDOUBLE_WORDS(jx,i0,i1,x);
00060        ix = jx&0x7fff;
00061 
00062     /* x is INF or NaN */
00063        if(ix==0x7fff) return x+x;
00064 
00065        h = 0.5;
00066        if (jx & 0x8000) h = -h;
00067     /* |x| in [0,25], return sign(x)*0.5*(E+E/(E+1))) */
00068        if (ix < 0x4003 || (ix == 0x4003 && i0 <= 0xc8000000)) { /* |x|<25 */
00069            if (ix<0x3fdf)          /* |x|<2**-32 */
00070               if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
00071            t = __expm1l(fabsl(x));
00072            if(ix<0x3fff) return h*(2.0*t-t*t/(t+one));
00073            return h*(t+t/(t+one));
00074        }
00075 
00076     /* |x| in [25, log(maxdouble)] return 0.5*exp(|x|) */
00077        if (ix < 0x400c || (ix == 0x400c && i0 < 0xb17217f7))
00078               return h*__ieee754_expl(fabsl(x));
00079 
00080     /* |x| in [log(maxdouble), overflowthreshold] */
00081        if (ix<0x400c || (ix == 0x400c && (i0 < 0xb174ddc0
00082                                       || (i0 == 0xb174ddc0
00083                                           && i1 <= 0x31aec0ea)))) {
00084            w = __ieee754_expl(0.5*fabsl(x));
00085            t = h*w;
00086            return t*w;
00087        }
00088 
00089     /* |x| > overflowthreshold, sinhl(x) overflow */
00090        return x*shuge;
00091 }