Line data    Source code

       1             : 
       2             : /* @(#)e_sinh.c 1.3 95/01/18 */
       3             : /*
       4             :  * ====================================================
       5             :  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
       6             :  *
       7             :  * Developed at SunSoft, a Sun Microsystems, Inc. business.
       8             :  * Permission to use, copy, modify, and distribute this
       9             :  * software is freely granted, provided that this notice 
      10             :  * is preserved.
      11             :  * ====================================================
      12             :  */
      13             : 
      14             : //#include <sys/cdefs.h>
      15             : //__FBSDID("$FreeBSD$");
      16             : 
      17             : /* __ieee754_sinh(x)
      18             :  * Method : 
      19             :  * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
      20             :  *      1. Replace x by |x| (sinh(-x) = -sinh(x)). 
      21             :  *      2. 
      22             :  *                                                  E + E/(E+1)
      23             :  *          0        <= x <= 22     :  sinh(x) := --------------, E=expm1(x)
      24             :  *                                                      2
      25             :  *
      26             :  *          22       <= x <= lnovft :  sinh(x) := exp(x)/2 
      27             :  *          lnovft   <= x <= ln2ovft:  sinh(x) := exp(x/2)/2 * exp(x/2)
      28             :  *          ln2ovft  <  x        :  sinh(x) := x*shuge (overflow)
      29             :  *
      30             :  * Special cases:
      31             :  *      sinh(x) is |x| if x is +INF, -INF, or NaN.
      32             :  *      only sinh(0)=0 is exact for finite x.
      33             :  */
      34             : 
      35             : #include <float.h>
      36             : 
      37             : #include "math_private.h"
      38             : 
      39             : static const double one = 1.0, shuge = 1.0e307;
      40             : 
      41             : double
      42           0 : __ieee754_sinh(double x)
      43             : {
      44             :         double t,h;
      45             :         int32_t ix,jx;
      46             : 
      47             :     /* High word of |x|. */
      48           0 :         GET_HIGH_WORD(jx,x);
      49           0 :         ix = jx&0x7fffffff;
      50             : 
      51             :     /* x is INF or NaN */
      52           0 :         if(ix>=0x7ff00000) return x+x;       
      53             : 
      54           0 :         h = 0.5;
      55           0 :         if (jx<0) h = -h;
      56             :     /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
      57           0 :         if (ix < 0x40360000) {               /* |x|<22 */
      58           0 :             if (ix<0x3e300000)               /* |x|<2**-28 */
      59           0 :                 if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
      60           0 :             t = expm1(fabs(x));
      61           0 :             if(ix<0x3ff00000) return h*(2.0*t-t*t/(t+one));
      62           0 :             return h*(t+t/(t+one));
      63             :         }
      64             : 
      65             :     /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
      66           0 :         if (ix < 0x40862E42)  return h*__ieee754_exp(fabs(x));
      67             : 
      68             :     /* |x| in [log(maxdouble), overflowthresold] */
      69           0 :         if (ix<=0x408633CE)
      70           0 :             return h*2.0*__ldexp_exp(fabs(x), -1);
      71             : 
      72             :     /* |x| > overflowthresold, sinh(x) overflow */
      73           0 :         return x*shuge;
      74             : }