nsz Git - libm/blob - src/math/coshl.c

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