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