1/* mpfr_tanh -- hyperbolic tangent
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
26int
27mpfr_tanh (mpfr_ptr y, mpfr_srcptr xt , mpfr_rnd_t rnd_mode)
28{
29  /****** Declaration ******/
30  mpfr_t x;
31  int inexact;
32  MPFR_SAVE_EXPO_DECL (expo);
33
34  MPFR_LOG_FUNC
35    (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (xt), mpfr_log_prec, xt, rnd_mode),
36     ("y[%Pu]=%.*Rg inexact=%d",
37      mpfr_get_prec (y), mpfr_log_prec, y, inexact));
38
39  /* Special value checking */
40  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt)))
41    {
42      if (MPFR_IS_NAN (xt))
43        {
44          MPFR_SET_NAN (y);
45          MPFR_RET_NAN;
46        }
47      else if (MPFR_IS_INF (xt))
48        {
49          /* tanh(inf) = 1 && tanh(-inf) = -1 */
50          return mpfr_set_si (y, MPFR_INT_SIGN (xt), rnd_mode);
51        }
52      else /* tanh (0) = 0 and xt is zero */
53        {
54          MPFR_ASSERTD (MPFR_IS_ZERO(xt));
55          MPFR_SET_ZERO (y);
56          MPFR_SET_SAME_SIGN (y, xt);
57          MPFR_RET (0);
58        }
59    }
60
61  /* tanh(x) = x - x^3/3 + ... so the error is < 2^(3*EXP(x)-1) */
62  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2 * MPFR_GET_EXP (xt), 1, 0,
63                                    rnd_mode, {});
64
65  MPFR_TMP_INIT_ABS (x, xt);
66
67  MPFR_SAVE_EXPO_MARK (expo);
68
69  /* General case */
70  {
71    /* Declaration of the intermediary variable */
72    mpfr_t t, te;
73    mpfr_exp_t d;
74
75    /* Declaration of the size variable */
76    mpfr_prec_t Ny = MPFR_PREC(y);   /* target precision */
77    mpfr_prec_t Nt;                  /* working precision */
78    long int err;                  /* error */
79    int sign = MPFR_SIGN (xt);
80    MPFR_ZIV_DECL (loop);
81    MPFR_GROUP_DECL (group);
82
83    /* First check for BIG overflow of exp(2*x):
84       For x > 0, exp(2*x) > 2^(2*x)
85       If 2 ^(2*x) > 2^emax or x>emax/2, there is an overflow */
86    if (MPFR_UNLIKELY (mpfr_cmp_si (x, __gmpfr_emax/2) >= 0)) {
87      /* initialise of intermediary variables
88         since 'set_one' label assumes the variables have been
89         initialize */
90      MPFR_GROUP_INIT_2 (group, MPFR_PREC_MIN, t, te);
91      goto set_one;
92    }
93
94    /* Compute the precision of intermediary variable */
95    /* The optimal number of bits: see algorithms.tex */
96    Nt = Ny + MPFR_INT_CEIL_LOG2 (Ny) + 4;
97    /* if x is small, there will be a cancellation in exp(2x)-1 */
98    if (MPFR_GET_EXP (x) < 0)
99      Nt += -MPFR_GET_EXP (x);
100
101    /* initialise of intermediary variable */
102    MPFR_GROUP_INIT_2 (group, Nt, t, te);
103
104    MPFR_ZIV_INIT (loop, Nt);
105    for (;;) {
106      /* tanh = (exp(2x)-1)/(exp(2x)+1) */
107      mpfr_mul_2ui (te, x, 1, MPFR_RNDN);  /* 2x */
108      /* since x > 0, we can only have an overflow */
109      mpfr_exp (te, te, MPFR_RNDN);        /* exp(2x) */
110      if (MPFR_UNLIKELY (MPFR_IS_INF (te))) {
111      set_one:
112        inexact = MPFR_FROM_SIGN_TO_INT (sign);
113        mpfr_set4 (y, __gmpfr_one, MPFR_RNDN, sign);
114        if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG_SIGN (sign)))
115          {
116            inexact = -inexact;
117            mpfr_nexttozero (y);
118          }
119        break;
120      }
121      d = MPFR_GET_EXP (te);              /* For Error calculation */
122      mpfr_add_ui (t, te, 1, MPFR_RNDD);   /* exp(2x) + 1*/
123      mpfr_sub_ui (te, te, 1, MPFR_RNDU);  /* exp(2x) - 1*/
124      d = d - MPFR_GET_EXP (te);
125      mpfr_div (t, te, t, MPFR_RNDN);      /* (exp(2x)-1)/(exp(2x)+1)*/
126
127      /* Calculation of the error */
128      d = MAX(3, d + 1);
129      err = Nt - (d + 1);
130
131      if (MPFR_LIKELY ((d <= Nt / 2) && MPFR_CAN_ROUND (t, err, Ny, rnd_mode)))
132        {
133          inexact = mpfr_set4 (y, t, rnd_mode, sign);
134          break;
135        }
136
137      /* if t=1, we still can round since |sinh(x)| < 1 */
138      if (MPFR_GET_EXP (t) == 1)
139        goto set_one;
140
141      /* Actualisation of the precision */
142      MPFR_ZIV_NEXT (loop, Nt);
143      MPFR_GROUP_REPREC_2 (group, Nt, t, te);
144    }
145    MPFR_ZIV_FREE (loop);
146    MPFR_GROUP_CLEAR (group);
147  }
148  MPFR_SAVE_EXPO_FREE (expo);
149  inexact = mpfr_check_range (y, inexact, rnd_mode);
150
151  return inexact;
152}
153
154