Branch data     Line data    Source code

       1                 :             : /* @(#)e_cosh.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                 :             : /* coshl(x)
      14                 :             :  * Method :
      15                 :             :  * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2
      16                 :             :  *      1. Replace x by |x| (coshl(x) = coshl(-x)).
      17                 :             :  *      2.
      18                 :             :  *                                                      [ exp(x) - 1 ]^2
      19                 :             :  *          0        <= x <= ln2/2  :  coshl(x) := 1 + -------------------
      20                 :             :  *                                                         2*exp(x)
      21                 :             :  *
      22                 :             :  *                                                 exp(x) +  1/exp(x)
      23                 :             :  *          ln2/2    <= x <= 22     :  coshl(x) := -------------------
      24                 :             :  *                                                         2
      25                 :             :  *          22       <= x <= lnovft :  coshl(x) := expl(x)/2
      26                 :             :  *          lnovft   <= x <= ln2ovft:  coshl(x) := expl(x/2)/2 * expl(x/2)
      27                 :             :  *          ln2ovft  <  x        :  coshl(x) := huge*huge (overflow)
      28                 :             :  *
      29                 :             :  * Special cases:
      30                 :             :  *      coshl(x) is |x| if x is +INF, -INF, or NaN.
      31                 :             :  *      only coshl(0)=1 is exact for finite x.
      32                 :             :  */
      33                 :             : 
      34                 :             : #include <openlibm_math.h>
      35                 :             : 
      36                 :             : #include "math_private.h"
      37                 :             : 
      38                 :             : static const long double one = 1.0, half=0.5, huge = 1.0e4900L;
      39                 :             : 
      40                 :             : long double
      41                 :           0 : coshl(long double x)
      42                 :             : {
      43                 :             :         long double t,w;
      44                 :             :         int32_t ex;
      45                 :             :         u_int32_t mx,lx;
      46                 :             : 
      47                 :             :     /* High word of |x|. */
      48                 :           0 :         GET_LDOUBLE_WORDS(ex,mx,lx,x);
      49                 :           0 :         ex &= 0x7fff;
      50                 :             : 
      51                 :             :     /* x is INF or NaN */
      52         [ #  # ]:           0 :         if(ex==0x7fff) return x*x;
      53                 :             : 
      54                 :             :     /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
      55   [ #  #  #  #  :           0 :         if(ex < 0x3ffd || (ex == 0x3ffd && mx < 0xb17217f7u)) {
                   #  # ]
      56                 :           0 :             t = expm1l(fabsl(x));
      57                 :           0 :             w = one+t;
      58         [ #  # ]:           0 :             if (ex<0x3fbc) return w; /* cosh(tiny) = 1 */
      59                 :           0 :             return one+(t*t)/(w+w);
      60                 :             :         }
      61                 :             : 
      62                 :             :     /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
      63   [ #  #  #  #  :           0 :         if (ex < 0x4003 || (ex == 0x4003 && mx < 0xb0000000u)) {
                   #  # ]
      64                 :           0 :                 t = expl(fabsl(x));
      65                 :           0 :                 return half*t+half/t;
      66                 :             :         }
      67                 :             : 
      68                 :             :     /* |x| in [22, ln(maxdouble)] return half*exp(|x|) */
      69   [ #  #  #  #  :           0 :         if (ex < 0x400c || (ex == 0x400c && mx < 0xb1700000u))
                   #  # ]
      70                 :           0 :                 return half*expl(fabsl(x));
      71                 :             : 
      72                 :             :     /* |x| in [log(maxdouble), log(2*maxdouble)) */
      73   [ #  #  #  # ]:           0 :         if (ex == 0x400c && (mx < 0xb174ddc0u
      74   [ #  #  #  # ]:           0 :                              || (mx == 0xb174ddc0u && lx < 0x31aec0ebu)))
      75                 :             :         {
      76                 :           0 :             w = expl(half*fabsl(x));
      77                 :           0 :             t = half*w;
      78                 :           0 :             return t*w;
      79                 :             :         }
      80                 :             : 
      81                 :             :     /* |x| >= log(2*maxdouble), cosh(x) overflow */
      82                 :           0 :         return huge*huge;
      83                 :             : }