Back to index

glibc  2.9
e_sinhl.c
Go to the documentation of this file.
00001 /* e_sinhl.c -- long double version of e_sinh.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 /* Changes for 128-bit long double are
00018    Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
00019    and are incorporated herein by permission of the author.  The author 
00020    reserves the right to distribute this material elsewhere under different
00021    copying permissions.  These modifications are distributed here under 
00022    the following terms:
00023 
00024     This library is free software; you can redistribute it and/or
00025     modify it under the terms of the GNU Lesser General Public
00026     License as published by the Free Software Foundation; either
00027     version 2.1 of the License, or (at your option) any later version.
00028 
00029     This library is distributed in the hope that it will be useful,
00030     but WITHOUT ANY WARRANTY; without even the implied warranty of
00031     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00032     Lesser General Public License for more details.
00033 
00034     You should have received a copy of the GNU Lesser General Public
00035     License along with this library; if not, write to the Free Software
00036     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00037 
00038 /* __ieee754_sinhl(x)
00039  * Method :
00040  * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
00041  *      1. Replace x by |x| (sinhl(-x) = -sinhl(x)).
00042  *      2.
00043  *                                                   E + E/(E+1)
00044  *          0        <= x <= 25     :  sinhl(x) := --------------, E=expm1l(x)
00045  *                                                       2
00046  *
00047  *          25       <= x <= lnovft :  sinhl(x) := expl(x)/2
00048  *          lnovft   <= x <= ln2ovft:  sinhl(x) := expl(x/2)/2 * expl(x/2)
00049  *          ln2ovft  <  x           :  sinhl(x) := x*shuge (overflow)
00050  *
00051  * Special cases:
00052  *      sinhl(x) is |x| if x is +INF, -INF, or NaN.
00053  *      only sinhl(0)=0 is exact for finite x.
00054  */
00055 
00056 #include "math.h"
00057 #include "math_private.h"
00058 
00059 #ifdef __STDC__
00060 static const long double one = 1.0, shuge = 1.0e4931L,
00061 ovf_thresh = 1.1357216553474703894801348310092223067821E4L;
00062 #else
00063 static long double one = 1.0, shuge = 1.0e4931L,
00064 ovf_thresh = 1.1357216553474703894801348310092223067821E4L;
00065 #endif
00066 
00067 #ifdef __STDC__
00068 long double
00069 __ieee754_sinhl (long double x)
00070 #else
00071 long double
00072 __ieee754_sinhl (x)
00073      long double x;
00074 #endif
00075 {
00076   long double t, w, h;
00077   u_int32_t jx, ix;
00078   ieee854_long_double_shape_type u;
00079 
00080   /* Words of |x|. */
00081   u.value = x;
00082   jx = u.parts32.w0;
00083   ix = jx & 0x7fffffff;
00084 
00085   /* x is INF or NaN */
00086   if (ix >= 0x7fff0000)
00087     return x + x;
00088 
00089   h = 0.5;
00090   if (jx & 0x80000000)
00091     h = -h;
00092 
00093   /* Absolute value of x.  */
00094   u.parts32.w0 = ix;
00095 
00096   /* |x| in [0,40], return sign(x)*0.5*(E+E/(E+1))) */
00097   if (ix <= 0x40044000)
00098     {
00099       if (ix < 0x3fc60000) /* |x| < 2^-57 */
00100        if (shuge + x > one)
00101          return x;          /* sinh(tiny) = tiny with inexact */
00102       t = __expm1l (u.value);
00103       if (ix < 0x3fff0000)
00104        return h * (2.0 * t - t * t / (t + one));
00105       return h * (t + t / (t + one));
00106     }
00107 
00108   /* |x| in [40, log(maxdouble)] return 0.5*exp(|x|) */
00109   if (ix <= 0x400c62e3) /* 11356.375 */
00110     return h * __ieee754_expl (u.value);
00111 
00112   /* |x| in [log(maxdouble), overflowthreshold]
00113      Overflow threshold is log(2 * maxdouble).  */
00114   if (u.value <= ovf_thresh)
00115     {
00116       w = __ieee754_expl (0.5 * u.value);
00117       t = h * w;
00118       return t * w;
00119     }
00120 
00121   /* |x| > overflowthreshold, sinhl(x) overflow */
00122   return x * shuge;
00123 }