Home | History | Annotate | Line # | Download | only in src 
      1 /* mpfr_sinh_cosh -- hyperbolic sine and cosine
      2 
      3 Copyright 2001-2023 Free Software Foundation, Inc.
      4 Contributed by the AriC and Caramba 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 3 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.LESSER.  If not, see
     20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
     21 51 Franklin St, Fifth Floor, Boston, 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, mpfr_rnd_t rnd_mode)
     33 {
     34   mpfr_t x;
     35   int inexact_sh, inexact_ch;
     36 
     37   MPFR_ASSERTN (sh != ch);
     38 
     39   MPFR_LOG_FUNC
     40     (("x[%Pd]=%.*Rg rnd=%d",
     41       mpfr_get_prec (xt), mpfr_log_prec, xt, rnd_mode),
     42      ("sh[%Pd]=%.*Rg ch[%Pd]=%.*Rg",
     43       mpfr_get_prec (sh), mpfr_log_prec, sh,
     44       mpfr_get_prec (ch), mpfr_log_prec, ch));
     45 
     46   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt)))
     47     {
     48       if (MPFR_IS_NAN (xt))
     49         {
     50           MPFR_SET_NAN (ch);
     51           MPFR_SET_NAN (sh);
     52           MPFR_RET_NAN;
     53         }
     54       else if (MPFR_IS_INF (xt))
     55         {
     56           MPFR_SET_INF (sh);
     57           MPFR_SET_SAME_SIGN (sh, xt);
     58           MPFR_SET_INF (ch);
     59           MPFR_SET_POS (ch);
     60           MPFR_RET (0);
     61         }
     62       else /* xt is zero */
     63         {
     64           MPFR_ASSERTD (MPFR_IS_ZERO (xt));
     65           MPFR_SET_ZERO (sh);                   /* sinh(0) = 0 */
     66           MPFR_SET_SAME_SIGN (sh, xt);
     67           inexact_sh = 0;
     68           inexact_ch = mpfr_set_ui (ch, 1, rnd_mode); /* cosh(0) = 1 */
     69           return INEX(inexact_sh,inexact_ch);
     70         }
     71     }
     72 
     73   /* Warning: if we use MPFR_FAST_COMPUTE_IF_SMALL_INPUT here, make sure
     74      that the code also works in case of overlap (see sin_cos.c) */
     75 
     76   MPFR_TMP_INIT_ABS (x, xt);
     77 
     78   {
     79     mpfr_t s, c, ti;
     80     mpfr_exp_t d;
     81     mpfr_prec_t N;    /* Precision of the intermediary variables */
     82     long int err;    /* Precision of error */
     83     MPFR_ZIV_DECL (loop);
     84     MPFR_SAVE_EXPO_DECL (expo);
     85     MPFR_GROUP_DECL (group);
     86 
     87     MPFR_SAVE_EXPO_MARK (expo);
     88 
     89     /* compute the precision of intermediary variable */
     90     N = MPFR_PREC (ch);
     91     N = MAX (N, MPFR_PREC (sh));
     92     /* the optimal number of bits : see algorithms.ps */
     93     N = N + MPFR_INT_CEIL_LOG2 (N) + 4;
     94 
     95     /* initialize of intermediary variables */
     96     MPFR_GROUP_INIT_3 (group, N, s, c, ti);
     97 
     98     /* First computation of sinh_cosh */
     99     MPFR_ZIV_INIT (loop, N);
    100     for (;;)
    101       {
    102         MPFR_BLOCK_DECL (flags);
    103 
    104         /* compute sinh_cosh */
    105         MPFR_BLOCK (flags, mpfr_exp (s, x, MPFR_RNDD));
    106         if (MPFR_OVERFLOW (flags))
    107           /* exp(x) does overflow */
    108           {
    109             /* since cosh(x) >= exp(x), cosh(x) overflows too */
    110             inexact_ch = mpfr_overflow (ch, rnd_mode, MPFR_SIGN_POS);
    111             /* sinh(x) may be representable */
    112             inexact_sh = mpfr_sinh (sh, xt, rnd_mode);
    113             MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, MPFR_FLAGS_OVERFLOW);
    114             break;
    115           }
    116         d = MPFR_GET_EXP (s);
    117         mpfr_ui_div (ti, 1, s, MPFR_RNDU);  /* 1/exp(x) */
    118         mpfr_add (c, s, ti, MPFR_RNDU);     /* exp(x) + 1/exp(x) */
    119         mpfr_sub (s, s, ti, MPFR_RNDN);     /* exp(x) - 1/exp(x) */
    120         mpfr_div_2ui (c, c, 1, MPFR_RNDN);  /* 1/2(exp(x) + 1/exp(x)) */
    121         mpfr_div_2ui (s, s, 1, MPFR_RNDN);  /* 1/2(exp(x) - 1/exp(x)) */
    122 
    123         /* it may be that s is zero (in fact, it can only occur when exp(x)=1,
    124            and thus ti=1 too) */
    125         if (MPFR_IS_ZERO (s))
    126           err = N; /* double the precision */
    127         else
    128           {
    129             /* calculation of the error */
    130             d = d - MPFR_GET_EXP (s) + 2;
    131             /* error estimate: err = N-(__gmpfr_ceil_log2(1+pow(2,d)));*/
    132             err = N - (MAX (d, 0) + 1);
    133             if (MPFR_LIKELY (MPFR_CAN_ROUND (s, err, MPFR_PREC (sh),
    134                                              rnd_mode) &&               \
    135                              MPFR_CAN_ROUND (c, err, MPFR_PREC (ch),
    136                                              rnd_mode)))
    137               {
    138                 inexact_sh = mpfr_set4 (sh, s, rnd_mode, MPFR_SIGN (xt));
    139                 inexact_ch = mpfr_set (ch, c, rnd_mode);
    140                 break;
    141               }
    142           }
    143         /* actualization of the precision */
    144         N += err;
    145         MPFR_ZIV_NEXT (loop, N);
    146         MPFR_GROUP_REPREC_3 (group, N, s, c, ti);
    147       }
    148     MPFR_ZIV_FREE (loop);
    149     MPFR_GROUP_CLEAR (group);
    150     MPFR_SAVE_EXPO_FREE (expo);
    151   }
    152 
    153   /* now, let's raise the flags if needed */
    154   inexact_sh = mpfr_check_range (sh, inexact_sh, rnd_mode);
    155   inexact_ch = mpfr_check_range (ch, inexact_ch, rnd_mode);
    156 
    157   return INEX(inexact_sh,inexact_ch);
    158 }
    159