1 /* mpfr_sinh -- hyperbolic sine
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 /* The computation of sinh is done by
27 sinh(x) = 1/2 [e^(x)-e^(-x)] */
30 mpfr_sinh (mpfr_ptr y
, mpfr_srcptr xt
, mpfr_rnd_t rnd_mode
)
36 (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (xt
), mpfr_log_prec
, xt
, rnd_mode
),
37 ("y[%Pu]=%.*Rg inexact=%d",
38 mpfr_get_prec (y
), mpfr_log_prec
, y
, inexact
));
40 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (xt
)))
47 else if (MPFR_IS_INF (xt
))
50 MPFR_SET_SAME_SIGN (y
, xt
);
55 MPFR_ASSERTD (MPFR_IS_ZERO (xt
));
56 MPFR_SET_ZERO (y
); /* sinh(0) = 0 */
57 MPFR_SET_SAME_SIGN (y
, xt
);
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,
66 MPFR_TMP_INIT_ABS (x
, xt
);
71 mpfr_prec_t Nt
; /* Precision of the intermediary variable */
72 long int err
; /* Precision of error */
74 MPFR_SAVE_EXPO_DECL (expo
);
75 MPFR_GROUP_DECL (group
);
77 MPFR_SAVE_EXPO_MARK (expo
);
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
);
87 /* initialise of intermediary variables */
88 MPFR_GROUP_INIT_2 (group
, Nt
, t
, ti
);
90 /* First computation of sinh */
91 MPFR_ZIV_INIT (loop
, Nt
);
94 MPFR_BLOCK_DECL (flags
);
97 MPFR_BLOCK (flags
, mpfr_exp (t
, x
, MPFR_RNDD
));
98 if (MPFR_OVERFLOW (flags
))
99 /* exp(x) does overflow */
101 /* sinh(x) = 2 * sinh(x/2) * cosh(x/2) */
102 mpfr_div_2ui (ti
, x
, 1, MPFR_RNDD
); /* exact */
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)
110 inexact
= mpfr_overflow (y
, rnd_mode
, MPFR_SIGN (xt
));
111 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo
, MPFR_FLAGS_OVERFLOW
);
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
);
119 /* multiplication below, error(t) <= 5 ulp(t) */
120 MPFR_BLOCK (flags
, mpfr_mul (t
, t
, ti
, MPFR_RNDD
));
121 if (MPFR_OVERFLOW (flags
))
123 inexact
= mpfr_overflow (y
, rnd_mode
, MPFR_SIGN (xt
));
124 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo
, MPFR_FLAGS_OVERFLOW
);
128 /* doubling below, exact */
129 MPFR_BLOCK (flags
, mpfr_mul_2ui (t
, t
, 1, MPFR_RNDN
));
130 if (MPFR_OVERFLOW (flags
))
132 inexact
= mpfr_overflow (y
, rnd_mode
, MPFR_SIGN (xt
));
133 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo
, MPFR_FLAGS_OVERFLOW
);
137 /* we have lost at most 3 bits of precision */
139 if (MPFR_LIKELY (MPFR_CAN_ROUND (t
, err
, MPFR_PREC (y
),
142 inexact
= mpfr_set4 (y
, t
, rnd_mode
, MPFR_SIGN (xt
));
145 err
= Nt
; /* double the precision */
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)) */
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 */
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
),
167 inexact
= mpfr_set4 (y
, t
, rnd_mode
, MPFR_SIGN (xt
));
173 /* actualisation of the precision */
175 MPFR_ZIV_NEXT (loop
, Nt
);
176 MPFR_GROUP_REPREC_2 (group
, Nt
, t
, ti
);
178 MPFR_ZIV_FREE (loop
);
179 MPFR_GROUP_CLEAR (group
);
180 MPFR_SAVE_EXPO_FREE (expo
);
183 return mpfr_check_range (y
, inexact
, rnd_mode
);