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 }