contrib/mpfr/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 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 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, mp_rnd_t rnd_mode)
  43 {
  44   int inexact;
  45   mp_prec_t p, q;
  46   mpfr_t tmp1, tmp2;
  47   mp_limb_t *tmp1p, *tmp2p;
  48   MPFR_SAVE_EXPO_DECL (expo);
  49   MPFR_ZIV_DECL (loop);
  50   MPFR_TMP_DECL(marker);
  51
  52   MPFR_LOG_FUNC (("a[%#R]=%R rnd=%d", a, a, rnd_mode),
  53                  ("r[%#R]=%R inexact=%d", r, r, inexact));
  54
  55   /* Special cases */
  56   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a)))
  57     {
  58       /* If a is NaN, the result is NaN */
  59       if (MPFR_IS_NAN (a))
  60         {
  61           MPFR_SET_NAN (r);
  62           MPFR_RET_NAN;
  63         }
  64       /* check for infinity before zero */
  65       else if (MPFR_IS_INF (a))
  66         {
  67           if (MPFR_IS_NEG (a))
  68             /* log(-Inf) = NaN */
  69             {
  70               MPFR_SET_NAN (r);
  71               MPFR_RET_NAN;
  72             }
  73           else /* log(+Inf) = +Inf */
  74             {
  75               MPFR_SET_INF (r);
  76               MPFR_SET_POS (r);
  77               MPFR_RET (0);
  78             }
  79         }
  80       else /* a is zero */
  81         {
  82           MPFR_ASSERTD (MPFR_IS_ZERO (a));
  83           MPFR_SET_INF (r);
  84           MPFR_SET_NEG (r);
  85           MPFR_RET (0); /* log(0) is an exact -infinity */
  86         }
  87     }
  88   /* If a is negative, the result is NaN */
  89   else if (MPFR_UNLIKELY (MPFR_IS_NEG (a)))
  90     {
  91       MPFR_SET_NAN (r);
  92       MPFR_RET_NAN;
  93     }
  94   /* If a is 1, the result is 0 */
  95   else if (MPFR_UNLIKELY (MPFR_GET_EXP (a) == 1 && mpfr_cmp_ui (a, 1) == 0))
  96     {
  97       MPFR_SET_ZERO (r);
  98       MPFR_SET_POS (r);
  99       MPFR_RET (0); /* only "normal" case where the result is exact */
 100     }
 101
 102   q = MPFR_PREC (r);
 103
 104   /* use initial precision about q+lg(q)+5 */
 105   p = q + 5 + 2 * MPFR_INT_CEIL_LOG2 (q);
 106   /* % ~(mp_prec_t)BITS_PER_MP_LIMB  ;
 107      m=q; while (m) { p++; m >>= 1; }  */
 108   /* if (MPFR_LIKELY(p % BITS_PER_MP_LIMB != 0))
 109       p += BITS_PER_MP_LIMB - (p%BITS_PER_MP_LIMB); */
 110
 111   MPFR_TMP_MARK(marker);
 112   MPFR_SAVE_EXPO_MARK (expo);
 113
 114   MPFR_ZIV_INIT (loop, p);
 115   for (;;)
 116     {
 117       mp_size_t size;
 118       long m;
 119       mp_exp_t cancel;
 120
 121       /* Calculus of m (depends on p) */
 122       m = (p + 1) / 2 - MPFR_GET_EXP (a) + 1;
 123
 124       /* All the mpfr_t needed have a precision of p */
 125       size = (p-1)/BITS_PER_MP_LIMB+1;
 126       MPFR_TMP_INIT (tmp1p, tmp1, p, size);
 127       MPFR_TMP_INIT (tmp2p, tmp2, p, size);
 128
 129       mpfr_mul_2si (tmp2, a, m, GMP_RNDN);    /* s=a*2^m,        err<=1 ulp  */
 130       mpfr_div (tmp1, __gmpfr_four, tmp2, GMP_RNDN);/* 4/s,      err<=2 ulps */
 131       mpfr_agm (tmp2, __gmpfr_one, tmp1, GMP_RNDN); /* AG(1,4/s),err<=3 ulps */
 132       mpfr_mul_2ui (tmp2, tmp2, 1, GMP_RNDN); /* 2*AG(1,4/s),    err<=3 ulps */
 133       mpfr_const_pi (tmp1, GMP_RNDN);         /* compute pi,     err<=1ulp   */
 134       mpfr_div (tmp2, tmp1, tmp2, GMP_RNDN);  /* pi/2*AG(1,4/s), err<=5ulps  */
 135       mpfr_const_log2 (tmp1, GMP_RNDN);      /* compute log(2),  err<=1ulp   */
 136       mpfr_mul_si (tmp1, tmp1, m, GMP_RNDN); /* compute m*log(2),err<=2ulps  */
 137       mpfr_sub (tmp1, tmp2, tmp1, GMP_RNDN); /* log(a),    err<=7ulps+cancel */
 138
 139       if (MPFR_LIKELY (MPFR_IS_PURE_FP (tmp1) && MPFR_IS_PURE_FP (tmp2)))
 140         {
 141           cancel = MPFR_GET_EXP (tmp2) - MPFR_GET_EXP (tmp1);
 142           MPFR_LOG_MSG (("canceled bits=%ld\n", (long) cancel));
 143           MPFR_LOG_VAR (tmp1);
 144           if (MPFR_UNLIKELY (cancel < 0))
 145             cancel = 0;
 146
 147           /* we have 7 ulps of error from the above roundings,
 148              4 ulps from the 4/s^2 second order term,
 149              plus the canceled bits */
 150           if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp1, p-cancel-4, q, rnd_mode)))
 151             break;
 152
 153           /* VL: I think it is better to have an increment that it isn't
 154              too low; in particular, the increment must be positive even
 155              if cancel = 0 (can this occur?). */
 156           p += cancel >= 8 ? cancel : 8;
 157         }
 158       else
 159         {
 160           /* TODO: find why this case can occur and what is best to do
 161              with it. */
 162           p += 32;
 163         }
 164
 165       MPFR_ZIV_NEXT (loop, p);
 166     }
 167   MPFR_ZIV_FREE (loop);
 168   inexact = mpfr_set (r, tmp1, rnd_mode);
 169   /* We clean */
 170   MPFR_TMP_FREE(marker);
 171
 172   MPFR_SAVE_EXPO_FREE (expo);
 173   return mpfr_check_range (r, inexact, rnd_mode);
 174 }