1 /* mpfr_sinh_cosh -- hyperbolic sine and cosine
3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 #define INEXPOS(y) ((y) == 0 ? 0 : (((y) > 0) ? 1 : 2))
27 #define INEX(y,z) (INEXPOS(y) | (INEXPOS(z) << 2))
29 /* The computations are done by
30 cosh(x) = 1/2 [e^(x)+e^(-x)]
31 sinh(x) = 1/2 [e^(x)-e^(-x)]
32 Adapted from mpfr_sinh.c */
35 mpfr_sinh_cosh (mpfr_ptr sh
, mpfr_ptr ch
, mpfr_srcptr xt
, mpfr_rnd_t rnd_mode
)
38 int inexact_sh
, inexact_ch
;
40 MPFR_ASSERTN (sh
!= ch
);
43 (("x[%Pu]=%.*Rg rnd=%d",
44 mpfr_get_prec (xt
), mpfr_log_prec
, xt
, rnd_mode
),
45 ("sh[%Pu]=%.*Rg ch[%Pu]=%.*Rg",
46 mpfr_get_prec (sh
), mpfr_log_prec
, sh
,
47 mpfr_get_prec (ch
), mpfr_log_prec
, ch
));
49 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt
)))
57 else if (MPFR_IS_INF (xt
))
60 MPFR_SET_SAME_SIGN (sh
, xt
);
67 MPFR_ASSERTD (MPFR_IS_ZERO (xt
));
68 MPFR_SET_ZERO (sh
); /* sinh(0) = 0 */
69 MPFR_SET_SAME_SIGN (sh
, xt
);
71 inexact_ch
= mpfr_set_ui (ch
, 1, rnd_mode
); /* cosh(0) = 1 */
72 return INEX(inexact_sh
,inexact_ch
);
76 /* Warning: if we use MPFR_FAST_COMPUTE_IF_SMALL_INPUT here, make sure
77 that the code also works in case of overlap (see sin_cos.c) */
79 MPFR_TMP_INIT_ABS (x
, xt
);
84 mpfr_prec_t N
; /* Precision of the intermediary variables */
85 long int err
; /* Precision of error */
87 MPFR_SAVE_EXPO_DECL (expo
);
88 MPFR_GROUP_DECL (group
);
90 MPFR_SAVE_EXPO_MARK (expo
);
92 /* compute the precision of intermediary variable */
94 N
= MAX (N
, MPFR_PREC (sh
));
95 /* the optimal number of bits : see algorithms.ps */
96 N
= N
+ MPFR_INT_CEIL_LOG2 (N
) + 4;
98 /* initialise of intermediary variables */
99 MPFR_GROUP_INIT_3 (group
, N
, s
, c
, ti
);
101 /* First computation of sinh_cosh */
102 MPFR_ZIV_INIT (loop
, N
);
105 MPFR_BLOCK_DECL (flags
);
107 /* compute sinh_cosh */
108 MPFR_BLOCK (flags
, mpfr_exp (s
, x
, MPFR_RNDD
));
109 if (MPFR_OVERFLOW (flags
))
110 /* exp(x) does overflow */
112 /* since cosh(x) >= exp(x), cosh(x) overflows too */
113 inexact_ch
= mpfr_overflow (ch
, rnd_mode
, MPFR_SIGN_POS
);
114 /* sinh(x) may be representable */
115 inexact_sh
= mpfr_sinh (sh
, xt
, rnd_mode
);
116 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo
, MPFR_FLAGS_OVERFLOW
);
119 d
= MPFR_GET_EXP (s
);
120 mpfr_ui_div (ti
, 1, s
, MPFR_RNDU
); /* 1/exp(x) */
121 mpfr_add (c
, s
, ti
, MPFR_RNDU
); /* exp(x) + 1/exp(x) */
122 mpfr_sub (s
, s
, ti
, MPFR_RNDN
); /* exp(x) - 1/exp(x) */
123 mpfr_div_2ui (c
, c
, 1, MPFR_RNDN
); /* 1/2(exp(x) + 1/exp(x)) */
124 mpfr_div_2ui (s
, s
, 1, MPFR_RNDN
); /* 1/2(exp(x) - 1/exp(x)) */
126 /* it may be that s is zero (in fact, it can only occur when exp(x)=1,
127 and thus ti=1 too) */
128 if (MPFR_IS_ZERO (s
))
129 err
= N
; /* double the precision */
132 /* calculation of the error */
133 d
= d
- MPFR_GET_EXP (s
) + 2;
134 /* error estimate: err = N-(__gmpfr_ceil_log2(1+pow(2,d)));*/
135 err
= N
- (MAX (d
, 0) + 1);
136 if (MPFR_LIKELY (MPFR_CAN_ROUND (s
, err
, MPFR_PREC (sh
),
138 MPFR_CAN_ROUND (c
, err
, MPFR_PREC (ch
),
141 inexact_sh
= mpfr_set4 (sh
, s
, rnd_mode
, MPFR_SIGN (xt
));
142 inexact_ch
= mpfr_set (ch
, c
, rnd_mode
);
146 /* actualisation of the precision */
148 MPFR_ZIV_NEXT (loop
, N
);
149 MPFR_GROUP_REPREC_3 (group
, N
, s
, c
, ti
);
151 MPFR_ZIV_FREE (loop
);
152 MPFR_GROUP_CLEAR (group
);
153 MPFR_SAVE_EXPO_FREE (expo
);
156 /* now, let's raise the flags if needed */
157 inexact_sh
= mpfr_check_range (sh
, inexact_sh
, rnd_mode
);
158 inexact_ch
= mpfr_check_range (ch
, inexact_ch
, rnd_mode
);
160 return INEX(inexact_sh
,inexact_ch
);