1 /* mpfr_acosh -- inverse hyperbolic cosine
3 Copyright 2001-2015 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 acosh is done by *
27 * acosh= ln(x + sqrt(x^2-1)) */
30 mpfr_acosh (mpfr_ptr y
, mpfr_srcptr x
, mpfr_rnd_t rnd_mode
)
32 MPFR_SAVE_EXPO_DECL (expo
);
37 ("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x
), mpfr_log_prec
, x
, rnd_mode
),
38 ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y
), mpfr_log_prec
, y
,
41 /* Deal with special cases */
42 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
44 /* Nan, or zero or -Inf */
45 if (MPFR_IS_INF (x
) && MPFR_IS_POS (x
))
51 else /* Nan, or zero or -Inf */
57 comp
= mpfr_cmp_ui (x
, 1);
58 if (MPFR_UNLIKELY (comp
< 0))
63 else if (MPFR_UNLIKELY (comp
== 0))
65 MPFR_SET_ZERO (y
); /* acosh(1) = 0 */
69 MPFR_SAVE_EXPO_MARK (expo
);
73 /* Declaration of the intermediary variables */
75 /* Declaration of the size variables */
76 mpfr_prec_t Ny
= MPFR_PREC(y
); /* Precision of output variable */
77 mpfr_prec_t Nt
; /* Precision of the intermediary variable */
78 mpfr_exp_t err
, exp_te
, d
; /* Precision of error */
81 /* compute the precision of intermediary variable */
82 /* the optimal number of bits : see algorithms.tex */
83 Nt
= Ny
+ 4 + MPFR_INT_CEIL_LOG2 (Ny
);
85 /* initialization of intermediary variables */
88 /* First computation of acosh */
89 MPFR_ZIV_INIT (loop
, Nt
);
92 MPFR_BLOCK_DECL (flags
);
95 MPFR_BLOCK (flags
, mpfr_mul (t
, x
, x
, MPFR_RNDD
)); /* x^2 */
96 if (MPFR_OVERFLOW (flags
))
101 /* As x is very large and the precision is not too large, we
102 assume that we obtain the same result by evaluating ln(2x).
103 We need to compute ln(x) + ln(2) as 2x can overflow. TODO:
104 write a proof and add an MPFR_ASSERTN. */
105 mpfr_log (t
, x
, MPFR_RNDN
); /* err(log) < 1/2 ulp(t) */
106 pln2
= Nt
- MPFR_PREC_MIN
< MPFR_GET_EXP (t
) ?
107 MPFR_PREC_MIN
: Nt
- MPFR_GET_EXP (t
);
108 mpfr_init2 (ln2
, pln2
);
109 mpfr_const_log2 (ln2
, MPFR_RNDN
); /* err(ln2) < 1/2 ulp(t) */
110 mpfr_add (t
, t
, ln2
, MPFR_RNDN
); /* err <= 3/2 ulp(t) */
116 exp_te
= MPFR_GET_EXP (t
);
117 mpfr_sub_ui (t
, t
, 1, MPFR_RNDD
); /* x^2-1 */
118 if (MPFR_UNLIKELY (MPFR_IS_ZERO (t
)))
120 /* This means that x is very close to 1: x = 1 + t with
121 t < 2^(-Nt). We have: acosh(x) = sqrt(2t) (1 - eps(t))
122 with 0 < eps(t) < t / 12. */
123 mpfr_sub_ui (t
, x
, 1, MPFR_RNDD
); /* t = x - 1 */
124 mpfr_mul_2ui (t
, t
, 1, MPFR_RNDN
); /* 2t */
125 mpfr_sqrt (t
, t
, MPFR_RNDN
); /* sqrt(2t) */
130 d
= exp_te
- MPFR_GET_EXP (t
);
131 mpfr_sqrt (t
, t
, MPFR_RNDN
); /* sqrt(x^2-1) */
132 mpfr_add (t
, t
, x
, MPFR_RNDN
); /* sqrt(x^2-1)+x */
133 mpfr_log (t
, t
, MPFR_RNDN
); /* ln(sqrt(x^2-1)+x) */
135 /* error estimate -- see algorithms.tex */
136 err
= 3 + MAX (1, d
) - MPFR_GET_EXP (t
);
137 /* error is bounded by 1/2 + 2^err <= 2^(max(0,1+err)) */
138 err
= MAX (0, 1 + err
);
142 if (MPFR_LIKELY (MPFR_CAN_ROUND (t
, Nt
- err
, Ny
, rnd_mode
)))
145 /* reactualisation of the precision */
146 MPFR_ZIV_NEXT (loop
, Nt
);
147 mpfr_set_prec (t
, Nt
);
149 MPFR_ZIV_FREE (loop
);
151 inexact
= mpfr_set (y
, t
, rnd_mode
);
156 MPFR_SAVE_EXPO_FREE (expo
);
157 return mpfr_check_range (y
, inexact
, rnd_mode
);