contrib/mpfr/sinh.c

   1 /* mpfr_sinh -- hyperbolic sine
   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  /* The computation of sinh is done by
  27     sinh(x) = 1/2 [e^(x)-e^(-x)]          */
  28
  29 int
  30 mpfr_sinh (mpfr_ptr y, mpfr_srcptr xt, mp_rnd_t rnd_mode)
  31 {
  32   mpfr_t x;
  33   int inexact;
  34
  35   MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", xt, xt, rnd_mode),
  36                  ("y[%#R]=%R inexact=%d", y, y, inexact));
  37
  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           MPFR_SET_INF (y);
  48           MPFR_SET_SAME_SIGN (y, xt);
  49           MPFR_RET (0);
  50         }
  51       else /* xt is zero */
  52         {
  53           MPFR_ASSERTD (MPFR_IS_ZERO (xt));
  54           MPFR_SET_ZERO (y);   /* sinh(0) = 0 */
  55           MPFR_SET_SAME_SIGN (y, xt);
  56           MPFR_RET (0);
  57         }
  58     }
  59
  60   /* sinh(x) = x + x^3/6 + ... so the error is < 2^(3*EXP(x)-2) */
  61   MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, xt, -2 * MPFR_GET_EXP(xt), 2, 1,
  62                                     rnd_mode, {});
  63
  64   MPFR_TMP_INIT_ABS (x, xt);
  65
  66   {
  67     mpfr_t t, ti;
  68     mp_exp_t d;
  69     mp_prec_t Nt;    /* Precision of the intermediary variable */
  70     long int err;    /* Precision of error */
  71     MPFR_ZIV_DECL (loop);
  72     MPFR_SAVE_EXPO_DECL (expo);
  73     MPFR_GROUP_DECL (group);
  74
  75     MPFR_SAVE_EXPO_MARK (expo);
  76
  77     /* compute the precision of intermediary variable */
  78     Nt = MAX (MPFR_PREC (x), MPFR_PREC (y));
  79     /* the optimal number of bits : see algorithms.ps */
  80     Nt = Nt + MPFR_INT_CEIL_LOG2 (Nt) + 4;
  81     /* If x is near 0, exp(x) - 1/exp(x) = 2*x+x^3/3+O(x^5) */
  82     if (MPFR_GET_EXP (x) < 0)
  83       Nt -= 2*MPFR_GET_EXP (x);
  84
  85     /* initialise of intermediary variables */
  86     MPFR_GROUP_INIT_2 (group, Nt, t, ti);
  87
  88     /* First computation of sinh */
  89     MPFR_ZIV_INIT (loop, Nt);
  90     for (;;)
  91       {
  92         MPFR_BLOCK_DECL (flags);
  93
  94         /* compute sinh */
  95         MPFR_BLOCK (flags, mpfr_exp (t, x, GMP_RNDD));
  96         if (MPFR_OVERFLOW (flags))
  97           /* exp(x) does overflow */
  98           {
  99             /* sinh(x) = 2 * sinh(x/2) * cosh(x/2) */
 100             mpfr_div_2ui (ti, x, 1, GMP_RNDD); /* exact */
 101
 102             /* t <- cosh(x/2): error(t) <= 1 ulp(t) */
 103             MPFR_BLOCK (flags, mpfr_cosh (t, ti, GMP_RNDD));
 104             if (MPFR_OVERFLOW (flags))
 105               /* when x>1 we have |sinh(x)| >= cosh(x/2), so sinh(x)
 106                  overflows too */
 107               {
 108                 inexact = mpfr_overflow (y, rnd_mode, MPFR_SIGN (xt));
 109                 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
 110                 break;
 111               }
 112
 113             /* ti <- sinh(x/2): , error(ti) <= 1 ulp(ti)
 114                cannot overflow because 0 < sinh(x) < cosh(x) when x > 0 */
 115             mpfr_sinh (ti, ti, GMP_RNDD);
 116
 117             /* multiplication below, error(t) <= 5 ulp(t) */
 118             MPFR_BLOCK (flags, mpfr_mul (t, t, ti, GMP_RNDD));
 119             if (MPFR_OVERFLOW (flags))
 120               {
 121                 inexact = mpfr_overflow (y, rnd_mode, MPFR_SIGN (xt));
 122                 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
 123                 break;
 124               }
 125
 126             /* doubling below, exact */
 127             MPFR_BLOCK (flags, mpfr_mul_2ui (t, t, 1, GMP_RNDN));
 128             if (MPFR_OVERFLOW (flags))
 129               {
 130                 inexact = mpfr_overflow (y, rnd_mode, MPFR_SIGN (xt));
 131                 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
 132                 break;
 133               }
 134
 135             /* we have lost at most 3 bits of precision */
 136             err = Nt - 3;
 137             if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, MPFR_PREC (y),
 138                                              rnd_mode)))
 139               {
 140                 inexact = mpfr_set4 (y, t, rnd_mode, MPFR_SIGN (xt));
 141                 break;
 142               }
 143             err = Nt; /* double the precision */
 144           }
 145         else
 146           {
 147             d = MPFR_GET_EXP (t);
 148             mpfr_ui_div (ti, 1, t, GMP_RNDU); /* 1/exp(x) */
 149             mpfr_sub (t, t, ti, GMP_RNDN);    /* exp(x) - 1/exp(x) */
 150             mpfr_div_2ui (t, t, 1, GMP_RNDN);  /* 1/2(exp(x) - 1/exp(x)) */
 151
 152             /* it may be that t is zero (in fact, it can only occur when te=1,
 153                and thus ti=1 too) */
 154             if (MPFR_IS_ZERO (t))
 155               err = Nt; /* double the precision */
 156             else
 157               {
 158                 /* calculation of the error */
 159                 d = d - MPFR_GET_EXP (t) + 2;
 160                 /* error estimate: err = Nt-(__gmpfr_ceil_log2(1+pow(2,d)));*/
 161                 err = Nt - (MAX (d, 0) + 1);
 162                 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, MPFR_PREC (y),
 163                                                  rnd_mode)))
 164                   {
 165                     inexact = mpfr_set4 (y, t, rnd_mode, MPFR_SIGN (xt));
 166                     break;
 167                   }
 168               }
 169           }
 170
 171         /* actualisation of the precision */
 172         Nt += err;
 173         MPFR_ZIV_NEXT (loop, Nt);
 174         MPFR_GROUP_REPREC_2 (group, Nt, t, ti);
 175       }
 176     MPFR_ZIV_FREE (loop);
 177     MPFR_GROUP_CLEAR (group);
 178     MPFR_SAVE_EXPO_FREE (expo);
 179   }
 180
 181   return mpfr_check_range (y, inexact, rnd_mode);
 182 }