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 contributed by
  13    Stephen L. Moshier <moshier@na-net.ornl.gov> */
  14
  15 /* __ieee754_coshl(x)
  16  * Method :
  17  * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2
  18  *      1. Replace x by |x| (coshl(x) = coshl(-x)).
  19  *      2.
  20  *                                                      [ exp(x) - 1 ]^2
  21  *          0        <= x <= ln2/2  :  coshl(x) := 1 + -------------------
  22  *                                                         2*exp(x)
  23  *
  24  *                                                 exp(x) +  1/exp(x)
  25  *          ln2/2    <= x <= 22     :  coshl(x) := -------------------
  26  *                                                         2
  27  *          22       <= x <= lnovft :  coshl(x) := expl(x)/2
  28  *          lnovft   <= x <= ln2ovft:  coshl(x) := expl(x/2)/2 * expl(x/2)
  29  *          ln2ovft  <  x           :  coshl(x) := huge*huge (overflow)
  30  *
  31  * Special cases:
  32  *      coshl(x) is |x| if x is +INF, -INF, or NaN.
  33  *      only coshl(0)=1 is exact for finite x.
  34  */
  35
  36 #include "math.h"
  37 #include "math_private.h"
  38
  39 #ifdef __STDC__
  40 static const long double one = 1.0, half = 0.5, huge = 1.0e4900L,
  41 ovf_thresh = 1.1357216553474703894801348310092223067821E4L;
  42 #else
  43 static long double one = 1.0, half = 0.5, huge = 1.0e4900L,
  44 ovf_thresh = 1.1357216553474703894801348310092223067821E4L;
  45 #endif
  46
  47 #ifdef __STDC__
  48 long double
  49 __ieee754_coshl (long double x)
  50 #else
  51 long double
  52 __ieee754_coshl (x)
  53      long double x;
  54 #endif
  55 {
  56   long double t, w;
  57   int32_t ex;
  58   u_int32_t mx, lx;
  59   ieee854_long_double_shape_type u;
  60
  61   u.value = x;
  62   ex = u.parts32.w0 & 0x7fffffff;
  63
  64   /* Absolute value of x.  */
  65   u.parts32.w0 = ex;
  66
  67   /* x is INF or NaN */
  68   if (ex >= 0x7fff0000)
  69     return x * x;
  70
  71   /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
  72   if (ex < 0x3ffd62e4) /* 0.3465728759765625 */
  73     {
  74       t = __expm1l (u.value);
  75       w = one + t;
  76       if (ex < 0x3fc60000) /* |x| < 2^-57 */
  77         return w;               /* cosh(tiny) = 1 */
  78
  79       return one + (t * t) / (w + w);
  80     }
  81
  82   /* |x| in [0.5*ln2,40], return (exp(|x|)+1/exp(|x|)/2; */
  83   if (ex < 0x40044000)
  84     {
  85       t = __ieee754_expl (u.value);
  86       return half * t + half / t;
  87     }
  88
  89   /* |x| in [22, ln(maxdouble)] return half*exp(|x|) */
  90   if (ex <= 0x400c62e3) /* 11356.375 */
  91     return half * __ieee754_expl (u.value);
  92
  93   /* |x| in [log(maxdouble), overflowthresold] */
  94   if (u.value <= ovf_thresh)
  95     {
  96       w = __ieee754_expl (half * u.value);
  97       t = half * w;
  98       return t * w;
  99     }
 100
 101   /* |x| > overflowthresold, cosh(x) overflow */
 102   return huge * huge;
 103 }