contrib/mpfr/src/log.c

   1 /* mpfr_log -- natural logarithm of a floating-point number
   2
   3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Caramel 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 3 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.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. */
  22
  23 #define MPFR_NEED_LONGLONG_H
  24 #include "mpfr-impl.h"
  25
  26 /* The computation of log(x) is done using the formula :
  27      if we want p bits of the result,
  28
  29                        pi
  30           log(x) ~ ------------  -   m log 2
  31                     2 AG(1,4/s)
  32
  33      where s = x 2^m > 2^(p/2)
  34
  35      More precisely, if F(x) = int(1/sqrt(1-(1-x^2)*sin(t)^2), t=0..PI/2),
  36      then for s>=1.26 we have log(s) < F(4/s) < log(s)*(1+4/s^2)
  37      from which we deduce pi/2/AG(1,4/s)*(1-4/s^2) < log(s) < pi/2/AG(1,4/s)
  38      so the relative error 4/s^2 is < 4/2^p i.e. 4 ulps.
  39 */
  40
  41 int
  42 mpfr_log (mpfr_ptr r, mpfr_srcptr a, mpfr_rnd_t rnd_mode)
  43 {
  44   int inexact;
  45   mpfr_prec_t p, q;
  46   mpfr_t tmp1, tmp2;
  47   MPFR_SAVE_EXPO_DECL (expo);
  48   MPFR_ZIV_DECL (loop);
  49   MPFR_GROUP_DECL(group);
  50
  51   MPFR_LOG_FUNC
  52     (("a[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (a), mpfr_log_prec, a, rnd_mode),
  53      ("r[%Pu]=%.*Rg inexact=%d", mpfr_get_prec (r), mpfr_log_prec, r,
  54       inexact));
  55
  56   /* Special cases */
  57   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a)))
  58     {
  59       /* If a is NaN, the result is NaN */
  60       if (MPFR_IS_NAN (a))
  61         {
  62           MPFR_SET_NAN (r);
  63           MPFR_RET_NAN;
  64         }
  65       /* check for infinity before zero */
  66       else if (MPFR_IS_INF (a))
  67         {
  68           if (MPFR_IS_NEG (a))
  69             /* log(-Inf) = NaN */
  70             {
  71               MPFR_SET_NAN (r);
  72               MPFR_RET_NAN;
  73             }
  74           else /* log(+Inf) = +Inf */
  75             {
  76               MPFR_SET_INF (r);
  77               MPFR_SET_POS (r);
  78               MPFR_RET (0);
  79             }
  80         }
  81       else /* a is zero */
  82         {
  83           MPFR_ASSERTD (MPFR_IS_ZERO (a));
  84           MPFR_SET_INF (r);
  85           MPFR_SET_NEG (r);
  86           mpfr_set_divby0 ();
  87           MPFR_RET (0); /* log(0) is an exact -infinity */
  88         }
  89     }
  90   /* If a is negative, the result is NaN */
  91   else if (MPFR_UNLIKELY (MPFR_IS_NEG (a)))
  92     {
  93       MPFR_SET_NAN (r);
  94       MPFR_RET_NAN;
  95     }
  96   /* If a is 1, the result is 0 */
  97   else if (MPFR_UNLIKELY (MPFR_GET_EXP (a) == 1 && mpfr_cmp_ui (a, 1) == 0))
  98     {
  99       MPFR_SET_ZERO (r);
 100       MPFR_SET_POS (r);
 101       MPFR_RET (0); /* only "normal" case where the result is exact */
 102     }
 103
 104   q = MPFR_PREC (r);
 105
 106   /* use initial precision about q+lg(q)+5 */
 107   p = q + 5 + 2 * MPFR_INT_CEIL_LOG2 (q);
 108   /* % ~(mpfr_prec_t)GMP_NUMB_BITS  ;
 109      m=q; while (m) { p++; m >>= 1; }  */
 110   /* if (MPFR_LIKELY(p % GMP_NUMB_BITS != 0))
 111       p += GMP_NUMB_BITS - (p%GMP_NUMB_BITS); */
 112
 113   MPFR_SAVE_EXPO_MARK (expo);
 114   MPFR_GROUP_INIT_2 (group, p, tmp1, tmp2);
 115
 116   MPFR_ZIV_INIT (loop, p);
 117   for (;;)
 118     {
 119       long m;
 120       mpfr_exp_t cancel;
 121
 122       /* Calculus of m (depends on p) */
 123       m = (p + 1) / 2 - MPFR_GET_EXP (a) + 1;
 124
 125       mpfr_mul_2si (tmp2, a, m, MPFR_RNDN);    /* s=a*2^m,        err<=1 ulp  */
 126       mpfr_div (tmp1, __gmpfr_four, tmp2, MPFR_RNDN);/* 4/s,      err<=2 ulps */
 127       mpfr_agm (tmp2, __gmpfr_one, tmp1, MPFR_RNDN); /* AG(1,4/s),err<=3 ulps */
 128       mpfr_mul_2ui (tmp2, tmp2, 1, MPFR_RNDN); /* 2*AG(1,4/s),    err<=3 ulps */
 129       mpfr_const_pi (tmp1, MPFR_RNDN);         /* compute pi,     err<=1ulp   */
 130       mpfr_div (tmp2, tmp1, tmp2, MPFR_RNDN);  /* pi/2*AG(1,4/s), err<=5ulps  */
 131       mpfr_const_log2 (tmp1, MPFR_RNDN);      /* compute log(2),  err<=1ulp   */
 132       mpfr_mul_si (tmp1, tmp1, m, MPFR_RNDN); /* compute m*log(2),err<=2ulps  */
 133       mpfr_sub (tmp1, tmp2, tmp1, MPFR_RNDN); /* log(a),    err<=7ulps+cancel */
 134
 135       if (MPFR_LIKELY (MPFR_IS_PURE_FP (tmp1) && MPFR_IS_PURE_FP (tmp2)))
 136         {
 137           cancel = MPFR_GET_EXP (tmp2) - MPFR_GET_EXP (tmp1);
 138           MPFR_LOG_MSG (("canceled bits=%ld\n", (long) cancel));
 139           MPFR_LOG_VAR (tmp1);
 140           if (MPFR_UNLIKELY (cancel < 0))
 141             cancel = 0;
 142
 143           /* we have 7 ulps of error from the above roundings,
 144              4 ulps from the 4/s^2 second order term,
 145              plus the canceled bits */
 146           if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp1, p-cancel-4, q, rnd_mode)))
 147             break;
 148
 149           /* VL: I think it is better to have an increment that it isn't
 150              too low; in particular, the increment must be positive even
 151              if cancel = 0 (can this occur?). */
 152           p += cancel >= 8 ? cancel : 8;
 153         }
 154       else
 155         {
 156           /* TODO: find why this case can occur and what is best to do
 157              with it. */
 158           p += 32;
 159         }
 160
 161       MPFR_ZIV_NEXT (loop, p);
 162       MPFR_GROUP_REPREC_2 (group, p, tmp1, tmp2);
 163     }
 164   MPFR_ZIV_FREE (loop);
 165   inexact = mpfr_set (r, tmp1, rnd_mode);
 166   /* We clean */
 167   MPFR_GROUP_CLEAR (group);
 168
 169   MPFR_SAVE_EXPO_FREE (expo);
 170   return mpfr_check_range (r, inexact, rnd_mode);
 171 }