contrib/mpfr/acosh.c

   1 /* mpfr_acosh -- inverse hyperbolic cosine
   2
   3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   7
   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 2.1 of the License, or (at your
  11 option) any later version.
  12
  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.
  17
  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.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 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
  29 int
  30 mpfr_acosh (mpfr_ptr y, mpfr_srcptr x , mp_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     mp_prec_t Ny = MPFR_PREC(y);   /* Precision of output variable */
  75     mp_prec_t Nt;                  /* Precision of the intermediary variable */
  76     mp_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, GMP_RNDD));  /* x^2 */
  94         if (MPFR_OVERFLOW (flags))
  95           {
  96             mpfr_t ln2;
  97             mp_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, GMP_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, GMP_RNDN);  /* err(ln2) < 1/2 ulp(t) */
 108             mpfr_add (t, t, ln2, GMP_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, GMP_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, GMP_RNDD);   /* t = x - 1 */
 122                 mpfr_mul_2ui (t, t, 1, GMP_RNDN);  /* 2t */
 123                 mpfr_sqrt (t, t, GMP_RNDN);        /* sqrt(2t) */
 124                 err = 1;
 125               }
 126             else
 127               {
 128                 d = exp_te - MPFR_GET_EXP (t);
 129                 mpfr_sqrt (t, t, GMP_RNDN);        /* sqrt(x^2-1) */
 130                 mpfr_add (t, t, x, GMP_RNDN);      /* sqrt(x^2-1)+x */
 131                 mpfr_log (t, t, GMP_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 }