1/* mpfr_acosh -- inverse hyperbolic cosine 2 3Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4Contributed by the Arenaire and Cacao projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 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 acosh is done by * 27 * acosh= ln(x + sqrt(x^2-1)) */ 28 29int 30mpfr_acosh (mpfr_ptr y, mpfr_srcptr x , mpfr_rnd_t rnd_mode) 31{ 32 MPFR_SAVE_EXPO_DECL (expo); 33 int inexact; 34 int comp; 35 36 MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd_mode), 37 ("y[%#R]=%R inexact=%d", y, y, inexact)); 38 39 /* Deal with special cases */ 40 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 41 { 42 /* Nan, or zero or -Inf */ 43 if (MPFR_IS_INF (x) && MPFR_IS_POS (x)) 44 { 45 MPFR_SET_INF (y); 46 MPFR_SET_POS (y); 47 MPFR_RET (0); 48 } 49 else /* Nan, or zero or -Inf */ 50 { 51 MPFR_SET_NAN (y); 52 MPFR_RET_NAN; 53 } 54 } 55 comp = mpfr_cmp_ui (x, 1); 56 if (MPFR_UNLIKELY (comp < 0)) 57 { 58 MPFR_SET_NAN (y); 59 MPFR_RET_NAN; 60 } 61 else if (MPFR_UNLIKELY (comp == 0)) 62 { 63 MPFR_SET_ZERO (y); /* acosh(1) = 0 */ 64 MPFR_SET_POS (y); 65 MPFR_RET (0); 66 } 67 MPFR_SAVE_EXPO_MARK (expo); 68 69 /* General case */ 70 { 71 /* Declaration of the intermediary variables */ 72 mpfr_t t; 73 /* Declaration of the size variables */ 74 mpfr_prec_t Ny = MPFR_PREC(y); /* Precision of output variable */ 75 mpfr_prec_t Nt; /* Precision of the intermediary variable */ 76 mpfr_exp_t err, exp_te, d; /* Precision of error */ 77 MPFR_ZIV_DECL (loop); 78 79 /* compute the precision of intermediary variable */ 80 /* the optimal number of bits : see algorithms.tex */ 81 Nt = Ny + 4 + MPFR_INT_CEIL_LOG2 (Ny); 82 83 /* initialization of intermediary variables */ 84 mpfr_init2 (t, Nt); 85 86 /* First computation of acosh */ 87 MPFR_ZIV_INIT (loop, Nt); 88 for (;;) 89 { 90 MPFR_BLOCK_DECL (flags); 91 92 /* compute acosh */ 93 MPFR_BLOCK (flags, mpfr_mul (t, x, x, MPFR_RNDD)); /* x^2 */ 94 if (MPFR_OVERFLOW (flags)) 95 { 96 mpfr_t ln2; 97 mpfr_prec_t pln2; 98 99 /* As x is very large and the precision is not too large, we 100 assume that we obtain the same result by evaluating ln(2x). 101 We need to compute ln(x) + ln(2) as 2x can overflow. TODO: 102 write a proof and add an MPFR_ASSERTN. */ 103 mpfr_log (t, x, MPFR_RNDN); /* err(log) < 1/2 ulp(t) */ 104 pln2 = Nt - MPFR_PREC_MIN < MPFR_GET_EXP (t) ? 105 MPFR_PREC_MIN : Nt - MPFR_GET_EXP (t); 106 mpfr_init2 (ln2, pln2); 107 mpfr_const_log2 (ln2, MPFR_RNDN); /* err(ln2) < 1/2 ulp(t) */ 108 mpfr_add (t, t, ln2, MPFR_RNDN); /* err <= 3/2 ulp(t) */ 109 mpfr_clear (ln2); 110 err = 1; 111 } 112 else 113 { 114 exp_te = MPFR_GET_EXP (t); 115 mpfr_sub_ui (t, t, 1, MPFR_RNDD); /* x^2-1 */ 116 if (MPFR_UNLIKELY (MPFR_IS_ZERO (t))) 117 { 118 /* This means that x is very close to 1: x = 1 + t with 119 t < 2^(-Nt). We have: acosh(x) = sqrt(2t) (1 - eps(t)) 120 with 0 < eps(t) < t / 12. */ 121 mpfr_sub_ui (t, x, 1, MPFR_RNDD); /* t = x - 1 */ 122 mpfr_mul_2ui (t, t, 1, MPFR_RNDN); /* 2t */ 123 mpfr_sqrt (t, t, MPFR_RNDN); /* sqrt(2t) */ 124 err = 1; 125 } 126 else 127 { 128 d = exp_te - MPFR_GET_EXP (t); 129 mpfr_sqrt (t, t, MPFR_RNDN); /* sqrt(x^2-1) */ 130 mpfr_add (t, t, x, MPFR_RNDN); /* sqrt(x^2-1)+x */ 131 mpfr_log (t, t, MPFR_RNDN); /* ln(sqrt(x^2-1)+x) */ 132 133 /* error estimate -- see algorithms.tex */ 134 err = 3 + MAX (1, d) - MPFR_GET_EXP (t); 135 /* error is bounded by 1/2 + 2^err <= 2^(max(0,1+err)) */ 136 err = MAX (0, 1 + err); 137 } 138 } 139 140 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, Nt - err, Ny, rnd_mode))) 141 break; 142 143 /* reactualisation of the precision */ 144 MPFR_ZIV_NEXT (loop, Nt); 145 mpfr_set_prec (t, Nt); 146 } 147 MPFR_ZIV_FREE (loop); 148 149 inexact = mpfr_set (y, t, rnd_mode); 150 151 mpfr_clear (t); 152 } 153 154 MPFR_SAVE_EXPO_FREE (expo); 155 return mpfr_check_range (y, inexact, rnd_mode); 156} 157