sysdeps/ieee754/ldbl-128/e_coshl.c

   1 /*
   2  * ====================================================
   3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   4  *
   5  * Developed at SunPro, a Sun Microsystems, Inc. business.
   6  * Permission to use, copy, modify, and distribute this
   7  * software is freely granted, provided that this notice
   8  * is preserved.
   9  * ====================================================
  10  */
  11
  12 /* Changes for 128-bit long double are
  13    Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
  14    and are incorporated herein by permission of the author.  The author
  15    reserves the right to distribute this material elsewhere under different
  16    copying permissions.  These modifications are distributed here under
  17    the following terms:
  18
  19     This library is free software; you can redistribute it and/or
  20     modify it under the terms of the GNU Lesser General Public
  21     License as published by the Free Software Foundation; either
  22     version 2.1 of the License, or (at your option) any later version.
  23
  24     This library is distributed in the hope that it will be useful,
  25     but WITHOUT ANY WARRANTY; without even the implied warranty of
  26     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  27     Lesser General Public License for more details.
  28
  29     You should have received a copy of the GNU Lesser General Public
  30     License along with this library; if not, write to the Free Software
  31     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
  32
  33 /* __ieee754_coshl(x)
  34  * Method :
  35  * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2
  36  *      1. Replace x by |x| (coshl(x) = coshl(-x)).
  37  *      2.
  38  *                                                      [ exp(x) - 1 ]^2
  39  *          0        <= x <= ln2/2  :  coshl(x) := 1 + -------------------
  40  *                                                         2*exp(x)
  41  *
  42  *                                                 exp(x) +  1/exp(x)
  43  *          ln2/2    <= x <= 22     :  coshl(x) := -------------------
  44  *                                                         2
  45  *          22       <= x <= lnovft :  coshl(x) := expl(x)/2
  46  *          lnovft   <= x <= ln2ovft:  coshl(x) := expl(x/2)/2 * expl(x/2)
  47  *          ln2ovft  <  x           :  coshl(x) := huge*huge (overflow)
  48  *
  49  * Special cases:
  50  *      coshl(x) is |x| if x is +INF, -INF, or NaN.
  51  *      only coshl(0)=1 is exact for finite x.
  52  */
  53
  54 #include "math.h"
  55 #include "math_private.h"
  56
  57 #ifdef __STDC__
  58 static const long double one = 1.0, half = 0.5, huge = 1.0e4900L,
  59 ovf_thresh = 1.1357216553474703894801348310092223067821E4L;
  60 #else
  61 static long double one = 1.0, half = 0.5, huge = 1.0e4900L,
  62 ovf_thresh = 1.1357216553474703894801348310092223067821E4L;
  63 #endif
  64
  65 #ifdef __STDC__
  66 long double
  67 __ieee754_coshl (long double x)
  68 #else
  69 long double
  70 __ieee754_coshl (x)
  71      long double x;
  72 #endif
  73 {
  74   long double t, w;
  75   int32_t ex;
  76   ieee854_long_double_shape_type u;
  77
  78   u.value = x;
  79   ex = u.parts32.w0 & 0x7fffffff;
  80
  81   /* Absolute value of x.  */
  82   u.parts32.w0 = ex;
  83
  84   /* x is INF or NaN */
  85   if (ex >= 0x7fff0000)
  86     return x * x;
  87
  88   /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
  89   if (ex < 0x3ffd62e4) /* 0.3465728759765625 */
  90     {
  91       t = __expm1l (u.value);
  92       w = one + t;
  93       if (ex < 0x3fb80000) /* |x| < 2^-116 */
  94         return w;               /* cosh(tiny) = 1 */
  95
  96       return one + (t * t) / (w + w);
  97     }
  98
  99   /* |x| in [0.5*ln2,40], return (exp(|x|)+1/exp(|x|)/2; */
 100   if (ex < 0x40044000)
 101     {
 102       t = __ieee754_expl (u.value);
 103       return half * t + half / t;
 104     }
 105
 106   /* |x| in [22, ln(maxdouble)] return half*exp(|x|) */
 107   if (ex <= 0x400c62e3) /* 11356.375 */
 108     return half * __ieee754_expl (u.value);
 109
 110   /* |x| in [log(maxdouble), overflowthresold] */
 111   if (u.value <= ovf_thresh)
 112     {
 113       w = __ieee754_expl (half * u.value);
 114       t = half * w;
 115       return t * w;
 116     }
 117
 118   /* |x| > overflowthresold, cosh(x) overflow */
 119   return huge * huge;
 120 }