Back to index

glibc  2.9
e_coshl.c
Go to the documentation of this file.
00001 /*
00002  * ====================================================
00003  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00004  *
00005  * Developed at SunPro, a Sun Microsystems, Inc. business.
00006  * Permission to use, copy, modify, and distribute this
00007  * software is freely granted, provided that this notice
00008  * is preserved.
00009  * ====================================================
00010  */
00011 
00012 /* Changes for 128-bit long double are
00013    Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
00014    and are incorporated herein by permission of the author.  The author 
00015    reserves the right to distribute this material elsewhere under different
00016    copying permissions.  These modifications are distributed here under 
00017    the following terms:
00018 
00019     This library is free software; you can redistribute it and/or
00020     modify it under the terms of the GNU Lesser General Public
00021     License as published by the Free Software Foundation; either
00022     version 2.1 of the License, or (at your option) any later version.
00023 
00024     This library is distributed in the hope that it will be useful,
00025     but WITHOUT ANY WARRANTY; without even the implied warranty of
00026     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00027     Lesser General Public License for more details.
00028 
00029     You should have received a copy of the GNU Lesser General Public
00030     License along with this library; if not, write to the Free Software
00031     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00032 
00033 /* __ieee754_coshl(x)
00034  * Method :
00035  * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2
00036  *      1. Replace x by |x| (coshl(x) = coshl(-x)).
00037  *      2.
00038  *                                                      [ exp(x) - 1 ]^2
00039  *          0        <= x <= ln2/2  :  coshl(x) := 1 + -------------------
00040  *                                                         2*exp(x)
00041  *
00042  *                                                 exp(x) +  1/exp(x)
00043  *          ln2/2    <= x <= 22     :  coshl(x) := -------------------
00044  *                                                         2
00045  *          22       <= x <= lnovft :  coshl(x) := expl(x)/2
00046  *          lnovft   <= x <= ln2ovft:  coshl(x) := expl(x/2)/2 * expl(x/2)
00047  *          ln2ovft  <  x           :  coshl(x) := huge*huge (overflow)
00048  *
00049  * Special cases:
00050  *      coshl(x) is |x| if x is +INF, -INF, or NaN.
00051  *      only coshl(0)=1 is exact for finite x.
00052  */
00053 
00054 #include "math.h"
00055 #include "math_private.h"
00056 
00057 #ifdef __STDC__
00058 static const long double one = 1.0, half = 0.5, huge = 1.0e4900L,
00059 ovf_thresh = 1.1357216553474703894801348310092223067821E4L;
00060 #else
00061 static long double one = 1.0, half = 0.5, huge = 1.0e4900L,
00062 ovf_thresh = 1.1357216553474703894801348310092223067821E4L;
00063 #endif
00064 
00065 #ifdef __STDC__
00066 long double
00067 __ieee754_coshl (long double x)
00068 #else
00069 long double
00070 __ieee754_coshl (x)
00071      long double x;
00072 #endif
00073 {
00074   long double t, w;
00075   int32_t ex;
00076   ieee854_long_double_shape_type u;
00077 
00078   u.value = x;
00079   ex = u.parts32.w0 & 0x7fffffff;
00080 
00081   /* Absolute value of x.  */
00082   u.parts32.w0 = ex;
00083 
00084   /* x is INF or NaN */
00085   if (ex >= 0x7fff0000)
00086     return x * x;
00087 
00088   /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */
00089   if (ex < 0x3ffd62e4) /* 0.3465728759765625 */
00090     {
00091       t = __expm1l (u.value);
00092       w = one + t;
00093       if (ex < 0x3fb80000) /* |x| < 2^-116 */
00094        return w;            /* cosh(tiny) = 1 */
00095 
00096       return one + (t * t) / (w + w);
00097     }
00098 
00099   /* |x| in [0.5*ln2,40], return (exp(|x|)+1/exp(|x|)/2; */
00100   if (ex < 0x40044000)
00101     {
00102       t = __ieee754_expl (u.value);
00103       return half * t + half / t;
00104     }
00105 
00106   /* |x| in [22, ln(maxdouble)] return half*exp(|x|) */
00107   if (ex <= 0x400c62e3) /* 11356.375 */
00108     return half * __ieee754_expl (u.value);
00109 
00110   /* |x| in [log(maxdouble), overflowthresold] */
00111   if (u.value <= ovf_thresh)
00112     {
00113       w = __ieee754_expl (half * u.value);
00114       t = half * w;
00115       return t * w;
00116     }
00117 
00118   /* |x| > overflowthresold, cosh(x) overflow */
00119   return huge * huge;
00120 }