1/* mpfr_asinh -- inverse hyperbolic sine 2 3Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel 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 asinh is done by * 27 * asinh = ln(x + sqrt(x^2 + 1)) */ 28 29int 30mpfr_asinh (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 31{ 32 int inexact; 33 int signx, neg; 34 mpfr_prec_t Ny, Nt; 35 mpfr_t t; /* auxiliary variables */ 36 mpfr_exp_t err; 37 MPFR_SAVE_EXPO_DECL (expo); 38 MPFR_ZIV_DECL (loop); 39 40 MPFR_LOG_FUNC ( 41 ("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode), 42 ("y[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (y), mpfr_log_prec, y, 43 inexact)); 44 45 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 46 { 47 if (MPFR_IS_NAN (x)) 48 { 49 MPFR_SET_NAN (y); 50 MPFR_RET_NAN; 51 } 52 else if (MPFR_IS_INF (x)) 53 { 54 MPFR_SET_INF (y); 55 MPFR_SET_SAME_SIGN (y, x); 56 MPFR_RET (0); 57 } 58 else /* x is necessarily 0 */ 59 { 60 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 61 MPFR_SET_ZERO (y); /* asinh(0) = 0 */ 62 MPFR_SET_SAME_SIGN (y, x); 63 MPFR_RET (0); 64 } 65 } 66 67 /* asinh(x) = x - x^3/6 + ... so the error is < 2^(3*EXP(x)-2) */ 68 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (y, x, -2 * MPFR_GET_EXP (x), 2, 0, 69 rnd_mode, {}); 70 71 Ny = MPFR_PREC (y); /* Precision of output variable */ 72 73 signx = MPFR_SIGN (x); 74 neg = MPFR_IS_NEG (x); 75 76 /* General case */ 77 78 /* compute the precision of intermediary variable */ 79 /* the optimal number of bits : see algorithms.tex */ 80 Nt = Ny + 4 + MPFR_INT_CEIL_LOG2 (Ny); 81 82 MPFR_SAVE_EXPO_MARK (expo); 83 84 /* initialize intermediary variables */ 85 mpfr_init2 (t, Nt); 86 87 /* First computation of asinh */ 88 MPFR_ZIV_INIT (loop, Nt); 89 for (;;) 90 { 91 /* compute asinh */ 92 mpfr_mul (t, x, x, MPFR_RNDD); /* x^2 */ 93 mpfr_add_ui (t, t, 1, MPFR_RNDD); /* x^2+1 */ 94 mpfr_sqrt (t, t, MPFR_RNDN); /* sqrt(x^2+1) */ 95 (neg ? mpfr_sub : mpfr_add) (t, t, x, MPFR_RNDN); /* sqrt(x^2+1)+x */ 96 mpfr_log (t, t, MPFR_RNDN); /* ln(sqrt(x^2+1)+x)*/ 97 98 if (MPFR_LIKELY (MPFR_IS_PURE_FP (t))) 99 { 100 /* error estimate -- see algorithms.tex */ 101 err = Nt - (MAX (4 - MPFR_GET_EXP (t), 0) + 1); 102 if (MPFR_LIKELY (MPFR_IS_ZERO (t) 103 || MPFR_CAN_ROUND (t, err, Ny, rnd_mode))) 104 break; 105 } 106 107 /* actualisation of the precision */ 108 MPFR_ZIV_NEXT (loop, Nt); 109 mpfr_set_prec (t, Nt); 110 } 111 MPFR_ZIV_FREE (loop); 112 113 inexact = mpfr_set4 (y, t, rnd_mode, signx); 114 115 mpfr_clear (t); 116 117 MPFR_SAVE_EXPO_FREE (expo); 118 return mpfr_check_range (y, inexact, rnd_mode); 119} 120