contrib/mpfr/sinh_cosh.c

   1 /* mpfr_sinh_cosh -- hyperbolic sine and cosine
   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 computations are done by
  27     cosh(x) = 1/2 [e^(x)+e^(-x)]
  28     sinh(x) = 1/2 [e^(x)-e^(-x)]
  29     Adapted from mpfr_sinh.c     */
  30
  31 int
  32 mpfr_sinh_cosh (mpfr_ptr sh, mpfr_ptr ch, mpfr_srcptr xt, mp_rnd_t rnd_mode)
  33 {
  34   mpfr_t x;
  35   int inexact, inexact_sh, inexact_ch;
  36
  37   MPFR_ASSERTN (sh != ch);
  38
  39   MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", xt, xt, rnd_mode),
  40                  ("sh[%#R]=%R ch[%#R]=%R inexact=%d", sh, sh, ch, ch, inexact));
  41
  42   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt)))
  43     {
  44       if (MPFR_IS_NAN (xt))
  45         {
  46           MPFR_SET_NAN (ch);
  47           MPFR_SET_NAN (sh);
  48           MPFR_RET_NAN;
  49         }
  50       else if (MPFR_IS_INF (xt))
  51         {
  52           MPFR_SET_INF (sh);
  53           MPFR_SET_SAME_SIGN (sh, xt);
  54           MPFR_SET_INF (ch);
  55           MPFR_SET_POS (ch);
  56           MPFR_RET (0);
  57         }
  58       else /* xt is zero */
  59         {
  60           MPFR_ASSERTD (MPFR_IS_ZERO (xt));
  61           MPFR_SET_ZERO (sh);                   /* sinh(0) = 0 */
  62           MPFR_SET_SAME_SIGN (sh, xt);
  63           return mpfr_set_ui (ch, 1, rnd_mode); /* cosh(0) = 1 */
  64         }
  65     }
  66
  67   /* Warning: if we use MPFR_FAST_COMPUTE_IF_SMALL_INPUT here, make sure
  68      that the code also works in case of overlap (see sin_cos.c) */
  69
  70   MPFR_TMP_INIT_ABS (x, xt);
  71
  72   {
  73     mpfr_t s, c, ti;
  74     mp_exp_t d;
  75     mp_prec_t N;    /* Precision of the intermediary variables */
  76     long int err;    /* Precision of error */
  77     MPFR_ZIV_DECL (loop);
  78     MPFR_SAVE_EXPO_DECL (expo);
  79     MPFR_GROUP_DECL (group);
  80
  81     MPFR_SAVE_EXPO_MARK (expo);
  82
  83     /* compute the precision of intermediary variable */
  84     N = MPFR_PREC (ch);
  85     N = MAX (N, MPFR_PREC (sh));
  86     N = MAX (N, MPFR_PREC (x));
  87     /* the optimal number of bits : see algorithms.ps */
  88     N = N + MPFR_INT_CEIL_LOG2 (N) + 4;
  89
  90     /* initialise of intermediary variables */
  91     MPFR_GROUP_INIT_3 (group, N, s, c, ti);
  92
  93     /* First computation of sinh_cosh */
  94     MPFR_ZIV_INIT (loop, N);
  95     for (;;)
  96       {
  97         MPFR_BLOCK_DECL (flags);
  98
  99         /* compute sinh_cosh */
 100         MPFR_BLOCK (flags, mpfr_exp (s, x, GMP_RNDD));
 101         if (MPFR_OVERFLOW (flags))
 102           /* exp(x) does overflow */
 103           {
 104             /* since cosh(x) >= exp(x), cosh(x) overflows too */
 105             inexact_ch = mpfr_overflow (ch, rnd_mode, MPFR_SIGN_POS);
 106             /* sinh(x) may be representable */
 107             inexact_sh = mpfr_sinh (sh, xt, rnd_mode);
 108             MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
 109             break;
 110           }
 111         d = MPFR_GET_EXP (s);
 112         mpfr_ui_div (ti, 1, s, GMP_RNDU);  /* 1/exp(x) */
 113         mpfr_add (c, s, ti, GMP_RNDU);     /* exp(x) + 1/exp(x) */
 114         mpfr_sub (s, s, ti, GMP_RNDN);     /* exp(x) - 1/exp(x) */
 115         mpfr_div_2ui (c, c, 1, GMP_RNDN);  /* 1/2(exp(x) + 1/exp(x)) */
 116         mpfr_div_2ui (s, s, 1, GMP_RNDN);  /* 1/2(exp(x) - 1/exp(x)) */
 117
 118         /* it may be that s is zero (in fact, it can only occur when exp(x)=1,
 119            and thus ti=1 too) */
 120         if (MPFR_IS_ZERO (s))
 121           err = N; /* double the precision */
 122         else
 123           {
 124             /* calculation of the error */
 125             d = d - MPFR_GET_EXP (s) + 2;
 126             /* error estimate: err = N-(__gmpfr_ceil_log2(1+pow(2,d)));*/
 127             err = N - (MAX (d, 0) + 1);
 128             if (MPFR_LIKELY (MPFR_CAN_ROUND (s, err, MPFR_PREC (sh),
 129                                              rnd_mode) &&               \
 130                              MPFR_CAN_ROUND (c, err, MPFR_PREC (ch),
 131                                              rnd_mode)))
 132               {
 133                 inexact_sh = mpfr_set4 (sh, s, rnd_mode, MPFR_SIGN (xt));
 134                 inexact_ch = mpfr_set (ch, c, rnd_mode);
 135                 break;
 136               }
 137           }
 138         /* actualisation of the precision */
 139         N += err;
 140         MPFR_ZIV_NEXT (loop, N);
 141         MPFR_GROUP_REPREC_3 (group, N, s, c, ti);
 142       }
 143     MPFR_ZIV_FREE (loop);
 144     MPFR_GROUP_CLEAR (group);
 145     MPFR_SAVE_EXPO_FREE (expo);
 146   }
 147
 148   /* now, let's raise the flags if needed */
 149   inexact = mpfr_check_range (sh, inexact_sh, rnd_mode);
 150   inexact = mpfr_check_range (ch, inexact_ch, rnd_mode) || inexact;
 151
 152   return inexact;
 153 }