Branch data     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 "cdefs-compat.h"
      15                 :             : //__FBSDID("$FreeBSD: src/lib/msun/src/e_sinh.c,v 1.11 2011/10/21 06:28:47 das Exp $");
      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                 :             : #include <openlibm_math.h>
      37                 :             : 
      38                 :             : #include "math_private.h"
      39                 :             : 
      40                 :             : static const double one = 1.0, shuge = 1.0e307;
      41                 :             : 
      42                 :             : OLM_DLLEXPORT double
      43                 :           7 : __ieee754_sinh(double x)
      44                 :             : {
      45                 :             :         double t,h;
      46                 :             :         int32_t ix,jx;
      47                 :             : 
      48                 :             :     /* High word of |x|. */
      49                 :           7 :         GET_HIGH_WORD(jx,x);
      50                 :           7 :         ix = jx&0x7fffffff;
      51                 :             : 
      52                 :             :     /* x is INF or NaN */
      53         [ +  + ]:           7 :         if(ix>=0x7ff00000) return x+x;       
      54                 :             : 
      55                 :           4 :         h = 0.5;
      56         [ +  + ]:           4 :         if (jx<0) h = -h;
      57                 :             :     /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
      58         [ +  - ]:           4 :         if (ix < 0x40360000) {               /* |x|<22 */
      59         [ +  + ]:           4 :             if (ix<0x3e300000)               /* |x|<2**-28 */
      60         [ +  - ]:           2 :                 if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
      61                 :           2 :             t = expm1(fabs(x));
      62         [ +  - ]:           2 :             if(ix<0x3ff00000) return h*(2.0*t-t*t/(t+one));
      63                 :           0 :             return h*(t+t/(t+one));
      64                 :             :         }
      65                 :             : 
      66                 :             :     /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
      67         [ #  # ]:           0 :         if (ix < 0x40862E42)  return h*__ieee754_exp(fabs(x));
      68                 :             : 
      69                 :             :     /* |x| in [log(maxdouble), overflowthresold] */
      70         [ #  # ]:           0 :         if (ix<=0x408633CE)
      71                 :           0 :             return h*2.0*__ldexp_exp(fabs(x), -1);
      72                 :             : 
      73                 :             :     /* |x| > overflowthresold, sinh(x) overflow */
      74                 :           0 :         return x*shuge;
      75                 :             : }
      76                 :             : 
      77                 :             : #if (LDBL_MANT_DIG == 53)
      78                 :             : openlibm_weak_reference(sinh, sinhl);
      79                 :             : #endif