1/* mpfr_sinh_cosh -- hyperbolic sine and cosine
2
3Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4Contributed by the AriC and Caramel projects, INRIA.
5
6This file is part of the GNU MPFR Library.
7
8The GNU MPFR Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MPFR Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23#define MPFR_NEED_LONGLONG_H
24#include "mpfr-impl.h"
25
26#define INEXPOS(y) ((y) == 0 ? 0 : (((y) > 0) ? 1 : 2))
27#define INEX(y,z) (INEXPOS(y) | (INEXPOS(z) << 2))
28
29 /* The computations are done by
30    cosh(x) = 1/2 [e^(x)+e^(-x)]
31    sinh(x) = 1/2 [e^(x)-e^(-x)]
32    Adapted from mpfr_sinh.c     */
33
34int
35mpfr_sinh_cosh (mpfr_ptr sh, mpfr_ptr ch, mpfr_srcptr xt, mpfr_rnd_t rnd_mode)
36{
37  mpfr_t x;
38  int inexact_sh, inexact_ch;
39
40  MPFR_ASSERTN (sh != ch);
41
42  MPFR_LOG_FUNC
43    (("x[%Pu]=%.*Rg rnd=%d",
44      mpfr_get_prec (xt), mpfr_log_prec, xt, rnd_mode),
45     ("sh[%Pu]=%.*Rg ch[%Pu]=%.*Rg",
46      mpfr_get_prec (sh), mpfr_log_prec, sh,
47      mpfr_get_prec (ch), mpfr_log_prec, ch));
48
49  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt)))
50    {
51      if (MPFR_IS_NAN (xt))
52        {
53          MPFR_SET_NAN (ch);
54          MPFR_SET_NAN (sh);
55          MPFR_RET_NAN;
56        }
57      else if (MPFR_IS_INF (xt))
58        {
59          MPFR_SET_INF (sh);
60          MPFR_SET_SAME_SIGN (sh, xt);
61          MPFR_SET_INF (ch);
62          MPFR_SET_POS (ch);
63          MPFR_RET (0);
64        }
65      else /* xt is zero */
66        {
67          MPFR_ASSERTD (MPFR_IS_ZERO (xt));
68          MPFR_SET_ZERO (sh);                   /* sinh(0) = 0 */
69          MPFR_SET_SAME_SIGN (sh, xt);
70          inexact_sh = 0;
71          inexact_ch = mpfr_set_ui (ch, 1, rnd_mode); /* cosh(0) = 1 */
72          return INEX(inexact_sh,inexact_ch);
73        }
74    }
75
76  /* Warning: if we use MPFR_FAST_COMPUTE_IF_SMALL_INPUT here, make sure
77     that the code also works in case of overlap (see sin_cos.c) */
78
79  MPFR_TMP_INIT_ABS (x, xt);
80
81  {
82    mpfr_t s, c, ti;
83    mpfr_exp_t d;
84    mpfr_prec_t N;    /* Precision of the intermediary variables */
85    long int err;    /* Precision of error */
86    MPFR_ZIV_DECL (loop);
87    MPFR_SAVE_EXPO_DECL (expo);
88    MPFR_GROUP_DECL (group);
89
90    MPFR_SAVE_EXPO_MARK (expo);
91
92    /* compute the precision of intermediary variable */
93    N = MPFR_PREC (ch);
94    N = MAX (N, MPFR_PREC (sh));
95    /* the optimal number of bits : see algorithms.ps */
96    N = N + MPFR_INT_CEIL_LOG2 (N) + 4;
97
98    /* initialise of intermediary variables */
99    MPFR_GROUP_INIT_3 (group, N, s, c, ti);
100
101    /* First computation of sinh_cosh */
102    MPFR_ZIV_INIT (loop, N);
103    for (;;)
104      {
105        MPFR_BLOCK_DECL (flags);
106
107        /* compute sinh_cosh */
108        MPFR_BLOCK (flags, mpfr_exp (s, x, MPFR_RNDD));
109        if (MPFR_OVERFLOW (flags))
110          /* exp(x) does overflow */
111          {
112            /* since cosh(x) >= exp(x), cosh(x) overflows too */
113            inexact_ch = mpfr_overflow (ch, rnd_mode, MPFR_SIGN_POS);
114            /* sinh(x) may be representable */
115            inexact_sh = mpfr_sinh (sh, xt, rnd_mode);
116            MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
117            break;
118          }
119        d = MPFR_GET_EXP (s);
120        mpfr_ui_div (ti, 1, s, MPFR_RNDU);  /* 1/exp(x) */
121        mpfr_add (c, s, ti, MPFR_RNDU);     /* exp(x) + 1/exp(x) */
122        mpfr_sub (s, s, ti, MPFR_RNDN);     /* exp(x) - 1/exp(x) */
123        mpfr_div_2ui (c, c, 1, MPFR_RNDN);  /* 1/2(exp(x) + 1/exp(x)) */
124        mpfr_div_2ui (s, s, 1, MPFR_RNDN);  /* 1/2(exp(x) - 1/exp(x)) */
125
126        /* it may be that s is zero (in fact, it can only occur when exp(x)=1,
127           and thus ti=1 too) */
128        if (MPFR_IS_ZERO (s))
129          err = N; /* double the precision */
130        else
131          {
132            /* calculation of the error */
133            d = d - MPFR_GET_EXP (s) + 2;
134            /* error estimate: err = N-(__gmpfr_ceil_log2(1+pow(2,d)));*/
135            err = N - (MAX (d, 0) + 1);
136            if (MPFR_LIKELY (MPFR_CAN_ROUND (s, err, MPFR_PREC (sh),
137                                             rnd_mode) &&               \
138                             MPFR_CAN_ROUND (c, err, MPFR_PREC (ch),
139                                             rnd_mode)))
140              {
141                inexact_sh = mpfr_set4 (sh, s, rnd_mode, MPFR_SIGN (xt));
142                inexact_ch = mpfr_set (ch, c, rnd_mode);
143                break;
144              }
145          }
146        /* actualisation of the precision */
147        N += err;
148        MPFR_ZIV_NEXT (loop, N);
149        MPFR_GROUP_REPREC_3 (group, N, s, c, ti);
150      }
151    MPFR_ZIV_FREE (loop);
152    MPFR_GROUP_CLEAR (group);
153    MPFR_SAVE_EXPO_FREE (expo);
154  }
155
156  /* now, let's raise the flags if needed */
157  inexact_sh = mpfr_check_range (sh, inexact_sh, rnd_mode);
158  inexact_ch = mpfr_check_range (ch, inexact_ch, rnd_mode);
159
160  return INEX(inexact_sh,inexact_ch);
161}
162