contrib/mpfr/tanh.c

   1 /* mpfr_tanh -- hyperbolic tangent
   2
   3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   7
   8 The GNU MPFR Library is free software; you can redistribute it and/or modify
   9 it under the terms of the GNU Lesser General Public License as published by
  10 the Free Software Foundation; either version 2.1 of the License, or (at your
  11 option) any later version.
  12
  13 The GNU MPFR Library is distributed in the hope that it will be useful, but
  14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MPFR Library; see the file COPYING.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 MA 02110-1301, USA. */
  22
  23 #define MPFR_NEED_LONGLONG_H
  24 #include "mpfr-impl.h"
  25
  26 int
  27 mpfr_tanh (mpfr_ptr y, mpfr_srcptr xt , mp_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 (("x[%#R]=%R rnd=%d", xt, xt, rnd_mode),
  35                  ("y[%#R]=%R inexact=%d", y, y, inexact));
  36
  37   /* Special value checking */
  38   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt)))
  39     {
  40       if (MPFR_IS_NAN (xt))
  41         {
  42           MPFR_SET_NAN (y);
  43           MPFR_RET_NAN;
  44         }
  45       else if (MPFR_IS_INF (xt))
  46         {
  47           /* tanh(inf) = 1 && tanh(-inf) = -1 */
  48           return mpfr_set_si (y, MPFR_INT_SIGN (xt), rnd_mode);
  49         }
  50       else /* tanh (0) = 0 and xt is zero */
  51         {
  52           MPFR_ASSERTD (MPFR_IS_ZERO(xt));
  53           MPFR_SET_ZERO (y);
  54           MPFR_SET_SAME_SIGN (y, xt);
  55           MPFR_RET (0);
  56         }
  57     }
  58
  59   /* tanh(x) = x - x^3/3 + ... so the error is < 2^(3*EXP(x)-1) */
  60   MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2 * MPFR_GET_EXP (xt), 1, 0,
  61                                     rnd_mode, {});
  62
  63   MPFR_TMP_INIT_ABS (x, xt);
  64
  65   MPFR_SAVE_EXPO_MARK (expo);
  66
  67   /* General case */
  68   {
  69     /* Declaration of the intermediary variable */
  70     mpfr_t t, te;
  71     mp_exp_t d;
  72
  73     /* Declaration of the size variable */
  74     mp_prec_t Ny = MPFR_PREC(y);   /* target precision */
  75     mp_prec_t Nt;                  /* working precision */
  76     long int err;                  /* error */
  77     int sign = MPFR_SIGN (xt);
  78     MPFR_ZIV_DECL (loop);
  79     MPFR_GROUP_DECL (group);
  80
  81     /* First check for BIG overflow of exp(2*x):
  82        For x > 0, exp(2*x) > 2^(2*x)
  83        If 2 ^(2*x) > 2^emax or x>emax/2, there is an overflow */
  84     if (MPFR_UNLIKELY (mpfr_cmp_si (x, __gmpfr_emax/2) >= 0)) {
  85       /* initialise of intermediary variables
  86          since 'set_one' label assumes the variables have been
  87          initialize */
  88       MPFR_GROUP_INIT_2 (group, MPFR_PREC_MIN, t, te);
  89       goto set_one;
  90     }
  91
  92     /* Compute the precision of intermediary variable */
  93     /* The optimal number of bits: see algorithms.tex */
  94     Nt = Ny + MPFR_INT_CEIL_LOG2 (Ny) + 4;
  95     /* if x is small, there will be a cancellation in exp(2x)-1 */
  96     if (MPFR_GET_EXP (x) < 0)
  97       Nt += -MPFR_GET_EXP (x);
  98
  99     /* initialise of intermediary variable */
 100     MPFR_GROUP_INIT_2 (group, Nt, t, te);
 101
 102     MPFR_ZIV_INIT (loop, Nt);
 103     for (;;) {
 104       /* tanh = (exp(2x)-1)/(exp(2x)+1) */
 105       mpfr_mul_2ui (te, x, 1, GMP_RNDN);  /* 2x */
 106       /* since x > 0, we can only have an overflow */
 107       mpfr_exp (te, te, GMP_RNDN);        /* exp(2x) */
 108       if (MPFR_UNLIKELY (MPFR_IS_INF (te))) {
 109       set_one:
 110         inexact = MPFR_FROM_SIGN_TO_INT (sign);
 111         mpfr_set4 (y, __gmpfr_one, GMP_RNDN, sign);
 112         if (MPFR_IS_LIKE_RNDZ (rnd_mode, MPFR_IS_NEG_SIGN (sign)))
 113           {
 114             inexact = -inexact;
 115             mpfr_nexttozero (y);
 116           }
 117         break;
 118       }
 119       d = MPFR_GET_EXP (te);              /* For Error calculation */
 120       mpfr_add_ui (t, te, 1, GMP_RNDD);   /* exp(2x) + 1*/
 121       mpfr_sub_ui (te, te, 1, GMP_RNDU);  /* exp(2x) - 1*/
 122       d = d - MPFR_GET_EXP (te);
 123       mpfr_div (t, te, t, GMP_RNDN);      /* (exp(2x)-1)/(exp(2x)+1)*/
 124
 125       /* Calculation of the error */
 126       d = MAX(3, d + 1);
 127       err = Nt - (d + 1);
 128
 129       if (MPFR_LIKELY ((d <= Nt / 2) && MPFR_CAN_ROUND (t, err, Ny, rnd_mode)))
 130         {
 131           inexact = mpfr_set4 (y, t, rnd_mode, sign);
 132           break;
 133         }
 134
 135       /* if t=1, we still can round since |sinh(x)| < 1 */
 136       if (MPFR_GET_EXP (t) == 1)
 137         goto set_one;
 138
 139       /* Actualisation of the precision */
 140       MPFR_ZIV_NEXT (loop, Nt);
 141       MPFR_GROUP_REPREC_2 (group, Nt, t, te);
 142     }
 143     MPFR_ZIV_FREE (loop);
 144     MPFR_GROUP_CLEAR (group);
 145   }
 146   MPFR_SAVE_EXPO_FREE (expo);
 147   inexact = mpfr_check_range (y, inexact, rnd_mode);
 148
 149   return inexact;
 150 }
 151