Branch data     Line data    Source code

       1                 :             : /* @(#)e_sinh.c 5.1 93/09/24 */
       2                 :             : /*
       3                 :             :  * ====================================================
       4                 :             :  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
       5                 :             :  *
       6                 :             :  * Developed at SunPro, a Sun Microsystems, Inc. business.
       7                 :             :  * Permission to use, copy, modify, and distribute this
       8                 :             :  * software is freely granted, provided that this notice
       9                 :             :  * is preserved.
      10                 :             :  * ====================================================
      11                 :             :  */
      12                 :             : 
      13                 :             : /* sinhl(x)
      14                 :             :  * Method :
      15                 :             :  * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
      16                 :             :  *      1. Replace x by |x| (sinhl(-x) = -sinhl(x)).
      17                 :             :  *      2.
      18                 :             :  *                                                   E + E/(E+1)
      19                 :             :  *          0        <= x <= 25     :  sinhl(x) := --------------, E=expm1l(x)
      20                 :             :  *                                                       2
      21                 :             :  *
      22                 :             :  *          25       <= x <= lnovft :  sinhl(x) := expl(x)/2
      23                 :             :  *          lnovft   <= x <= ln2ovft:  sinhl(x) := expl(x/2)/2 * expl(x/2)
      24                 :             :  *          ln2ovft  <  x        :  sinhl(x) := x*shuge (overflow)
      25                 :             :  *
      26                 :             :  * Special cases:
      27                 :             :  *      sinhl(x) is |x| if x is +INF, -INF, or NaN.
      28                 :             :  *      only sinhl(0)=0 is exact for finite x.
      29                 :             :  */
      30                 :             : 
      31                 :             : #include <openlibm_math.h>
      32                 :             : 
      33                 :             : #include "math_private.h"
      34                 :             : 
      35                 :             : static const long double one = 1.0, shuge = 1.0e4931L;
      36                 :             : 
      37                 :             : long double
      38                 :           0 : sinhl(long double x)
      39                 :             : {
      40                 :             :         long double t,w,h;
      41                 :             :         u_int32_t jx,ix,i0,i1;
      42                 :             : 
      43                 :             :     /* Words of |x|. */
      44                 :           0 :         GET_LDOUBLE_WORDS(jx,i0,i1,x);
      45                 :           0 :         ix = jx&0x7fff;
      46                 :             : 
      47                 :             :     /* x is INF or NaN */
      48         [ #  # ]:           0 :         if(ix==0x7fff) return x+x;
      49                 :             : 
      50                 :           0 :         h = 0.5;
      51         [ #  # ]:           0 :         if (jx & 0x8000) h = -h;
      52                 :             :     /* |x| in [0,25], return sign(x)*0.5*(E+E/(E+1))) */
      53   [ #  #  #  #  :           0 :         if (ix < 0x4003 || (ix == 0x4003 && i0 <= 0xc8000000)) { /* |x|<25 */
                   #  # ]
      54         [ #  # ]:           0 :             if (ix<0x3fdf)            /* |x|<2**-32 */
      55         [ #  # ]:           0 :                 if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */
      56                 :           0 :             t = expm1l(fabsl(x));
      57         [ #  # ]:           0 :             if(ix<0x3fff) return h*(2.0*t-t*t/(t+one));
      58                 :           0 :             return h*(t+t/(t+one));
      59                 :             :         }
      60                 :             : 
      61                 :             :     /* |x| in [25, log(maxdouble)] return 0.5*exp(|x|) */
      62   [ #  #  #  #  :           0 :         if (ix < 0x400c || (ix == 0x400c && i0 < 0xb17217f7))
                   #  # ]
      63                 :           0 :                 return h*expl(fabsl(x));
      64                 :             : 
      65                 :             :     /* |x| in [log(maxdouble), overflowthreshold] */
      66   [ #  #  #  #  :           0 :         if (ix<0x400c || (ix == 0x400c && (i0 < 0xb174ddc0
                   #  # ]
      67         [ #  # ]:           0 :                                            || (i0 == 0xb174ddc0
      68         [ #  # ]:           0 :                                                && i1 <= 0x31aec0ea)))) {
      69                 :           0 :             w = expl(0.5*fabsl(x));
      70                 :           0 :             t = h*w;
      71                 :           0 :             return t*w;
      72                 :             :         }
      73                 :             : 
      74                 :             :     /* |x| > overflowthreshold, sinhl(x) overflow */
      75                 :           0 :         return x*shuge;
      76                 :             : }