152419Sjulian/* mpfr_acosh -- inverse hyperbolic cosine 252419Sjulian 352419SjulianCopyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 452419SjulianContributed by the AriC and Caramel projects, INRIA. 552419Sjulian 652419SjulianThis file is part of the GNU MPFR Library. 752419Sjulian 852419SjulianThe GNU MPFR Library is free software; you can redistribute it and/or modify 952419Sjulianit under the terms of the GNU Lesser General Public License as published by 1052419Sjulianthe Free Software Foundation; either version 3 of the License, or (at your 1152419Sjulianoption) any later version. 1252419Sjulian 1352419SjulianThe GNU MPFR Library is distributed in the hope that it will be useful, but 1452419SjulianWITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 1552419Sjulianor FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 1652419SjulianLicense for more details. 1752419Sjulian 1852419SjulianYou should have received a copy of the GNU Lesser General Public License 1952419Sjulianalong with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 2052419Sjulianhttp://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2152419Sjulian51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 2252419Sjulian 2352419Sjulian#define MPFR_NEED_LONGLONG_H 2452419Sjulian#include "mpfr-impl.h" 2552419Sjulian 2652419Sjulian/* The computation of acosh is done by * 2752419Sjulian * acosh= ln(x + sqrt(x^2-1)) */ 2852419Sjulian 2952419Sjulianint 3052419Sjulianmpfr_acosh (mpfr_ptr y, mpfr_srcptr x , mpfr_rnd_t rnd_mode) 3152419Sjulian{ 3252419Sjulian MPFR_SAVE_EXPO_DECL (expo); 3352419Sjulian int inexact; 3452419Sjulian int comp; 3552419Sjulian 3652419Sjulian MPFR_LOG_FUNC ( 3767506Sjulian ("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode), 3852419Sjulian ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, 3952419Sjulian inexact)); 4052752Sjulian 4152419Sjulian /* Deal with special cases */ 4252419Sjulian if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 4352419Sjulian { 4452419Sjulian /* Nan, or zero or -Inf */ 4552419Sjulian if (MPFR_IS_INF (x) && MPFR_IS_POS (x)) 4652419Sjulian { 4752419Sjulian MPFR_SET_INF (y); 4852419Sjulian MPFR_SET_POS (y); 4952419Sjulian MPFR_RET (0); 5052419Sjulian } 5152419Sjulian else /* Nan, or zero or -Inf */ 5252419Sjulian { 5352419Sjulian MPFR_SET_NAN (y); 5452419Sjulian MPFR_RET_NAN; 5552419Sjulian } 5652419Sjulian } 5752419Sjulian comp = mpfr_cmp_ui (x, 1); 5852419Sjulian if (MPFR_UNLIKELY (comp < 0)) 5952419Sjulian { 6052419Sjulian MPFR_SET_NAN (y); 6152419Sjulian MPFR_RET_NAN; 6252419Sjulian } 6352419Sjulian else if (MPFR_UNLIKELY (comp == 0)) 6452419Sjulian { 6552419Sjulian MPFR_SET_ZERO (y); /* acosh(1) = 0 */ 6652419Sjulian MPFR_SET_POS (y); 6752419Sjulian MPFR_RET (0); 6852419Sjulian } 6952919Sjulian MPFR_SAVE_EXPO_MARK (expo); 7052419Sjulian 7152419Sjulian /* General case */ 7252419Sjulian { 7352419Sjulian /* Declaration of the intermediary variables */ 7452419Sjulian mpfr_t t; 7552419Sjulian /* Declaration of the size variables */ 7652419Sjulian mpfr_prec_t Ny = MPFR_PREC(y); /* Precision of output variable */ 7752419Sjulian mpfr_prec_t Nt; /* Precision of the intermediary variable */ 7852419Sjulian mpfr_exp_t err, exp_te, d; /* Precision of error */ 7952419Sjulian MPFR_ZIV_DECL (loop); 8052419Sjulian 8152419Sjulian /* compute the precision of intermediary variable */ 8252419Sjulian /* the optimal number of bits : see algorithms.tex */ 8352419Sjulian Nt = Ny + 4 + MPFR_INT_CEIL_LOG2 (Ny); 8452419Sjulian 8552419Sjulian /* initialization of intermediary variables */ 8652419Sjulian mpfr_init2 (t, Nt); 8752419Sjulian 8852419Sjulian /* First computation of acosh */ 8952419Sjulian MPFR_ZIV_INIT (loop, Nt); 9052419Sjulian for (;;) 9152419Sjulian { 9252419Sjulian MPFR_BLOCK_DECL (flags); 9352419Sjulian 9452419Sjulian /* compute acosh */ 9552419Sjulian MPFR_BLOCK (flags, mpfr_mul (t, x, x, MPFR_RNDD)); /* x^2 */ 9652419Sjulian if (MPFR_OVERFLOW (flags)) 9752419Sjulian { 9852419Sjulian mpfr_t ln2; 9952752Sjulian mpfr_prec_t pln2; 10052752Sjulian 10152752Sjulian /* As x is very large and the precision is not too large, we 10252752Sjulian assume that we obtain the same result by evaluating ln(2x). 10352752Sjulian We need to compute ln(x) + ln(2) as 2x can overflow. TODO: 10452885Sjulian write a proof and add an MPFR_ASSERTN. */ 10552419Sjulian mpfr_log (t, x, MPFR_RNDN); /* err(log) < 1/2 ulp(t) */ 10652419Sjulian pln2 = Nt - MPFR_PREC_MIN < MPFR_GET_EXP (t) ? 10752419Sjulian MPFR_PREC_MIN : Nt - MPFR_GET_EXP (t); 10852419Sjulian mpfr_init2 (ln2, pln2); 10952419Sjulian mpfr_const_log2 (ln2, MPFR_RNDN); /* err(ln2) < 1/2 ulp(t) */ 11052419Sjulian mpfr_add (t, t, ln2, MPFR_RNDN); /* err <= 3/2 ulp(t) */ 11152419Sjulian mpfr_clear (ln2); 11252419Sjulian err = 1; 11352419Sjulian } 11452419Sjulian else 11552419Sjulian { 11652419Sjulian exp_te = MPFR_GET_EXP (t); 11752419Sjulian mpfr_sub_ui (t, t, 1, MPFR_RNDD); /* x^2-1 */ 11852419Sjulian if (MPFR_UNLIKELY (MPFR_IS_ZERO (t))) 11952419Sjulian { 12052419Sjulian /* This means that x is very close to 1: x = 1 + t with 12152419Sjulian t < 2^(-Nt). We have: acosh(x) = sqrt(2t) (1 - eps(t)) 12252419Sjulian with 0 < eps(t) < t / 12. */ 12370159Sjulian mpfr_sub_ui (t, x, 1, MPFR_RNDD); /* t = x - 1 */ 12452419Sjulian mpfr_mul_2ui (t, t, 1, MPFR_RNDN); /* 2t */ 12552419Sjulian mpfr_sqrt (t, t, MPFR_RNDN); /* sqrt(2t) */ 12652419Sjulian err = 1; 12752419Sjulian } 12852419Sjulian else 12952419Sjulian { 13052419Sjulian d = exp_te - MPFR_GET_EXP (t); 13152419Sjulian mpfr_sqrt (t, t, MPFR_RNDN); /* sqrt(x^2-1) */ 13252419Sjulian mpfr_add (t, t, x, MPFR_RNDN); /* sqrt(x^2-1)+x */ 13353913Sarchie mpfr_log (t, t, MPFR_RNDN); /* ln(sqrt(x^2-1)+x) */ 13453913Sarchie 13552419Sjulian /* error estimate -- see algorithms.tex */ 13652419Sjulian err = 3 + MAX (1, d) - MPFR_GET_EXP (t); 13752419Sjulian /* error is bounded by 1/2 + 2^err <= 2^(max(0,1+err)) */ 13852419Sjulian err = MAX (0, 1 + err); 13964512Sarchie } 14052419Sjulian } 14152419Sjulian 14252419Sjulian if (MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - err, Ny, rnd_mode))) 14360938Sjake break; 14452419Sjulian 14553526Sjulian /* reactualisation of the precision */ 14652419Sjulian MPFR_ZIV_NEXT (loop, Nt); 14752419Sjulian mpfr_set_prec (t, Nt); 14852419Sjulian } 14952419Sjulian MPFR_ZIV_FREE (loop); 15052419Sjulian 15152419Sjulian inexact = mpfr_set (y, t, rnd_mode); 15252419Sjulian 15352419Sjulian mpfr_clear (t); 15452419Sjulian } 15552419Sjulian 15652419Sjulian MPFR_SAVE_EXPO_FREE (expo); 15752419Sjulian return mpfr_check_range (y, inexact, rnd_mode); 15852419Sjulian} 15952419Sjulian