Home | History | Annotate | Line # | Download | only in src 
      1 /* mpfr_log -- natural logarithm of a floating-point number
      2 
      3 Copyright 1999-2023 Free Software Foundation, Inc.
      4 Contributed by the AriC and Caramba 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 https://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_exp_t exp_a;
     48   MPFR_SAVE_EXPO_DECL (expo);
     49   MPFR_ZIV_DECL (loop);
     50   MPFR_GROUP_DECL(group);
     51 
     52   MPFR_LOG_FUNC
     53     (("a[%Pd]=%.*Rg rnd=%d", mpfr_get_prec (a), mpfr_log_prec, a, rnd_mode),
     54      ("r[%Pd]=%.*Rg inexact=%d", mpfr_get_prec (r), mpfr_log_prec, r,
     55       inexact));
     56 
     57   /* Special cases */
     58   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a)))
     59     {
     60       /* If a is NaN, the result is NaN */
     61       if (MPFR_IS_NAN (a))
     62         {
     63           MPFR_SET_NAN (r);
     64           MPFR_RET_NAN;
     65         }
     66       /* check for infinity before zero */
     67       else if (MPFR_IS_INF (a))
     68         {
     69           if (MPFR_IS_NEG (a))
     70             /* log(-Inf) = NaN */
     71             {
     72               MPFR_SET_NAN (r);
     73               MPFR_RET_NAN;
     74             }
     75           else /* log(+Inf) = +Inf */
     76             {
     77               MPFR_SET_INF (r);
     78               MPFR_SET_POS (r);
     79               MPFR_RET (0);
     80             }
     81         }
     82       else /* a is zero */
     83         {
     84           MPFR_ASSERTD (MPFR_IS_ZERO (a));
     85           MPFR_SET_INF (r);
     86           MPFR_SET_NEG (r);
     87           MPFR_SET_DIVBY0 ();
     88           MPFR_RET (0); /* log(0) is an exact -infinity */
     89         }
     90     }
     91 
     92   /* If a is negative, the result is NaN */
     93   if (MPFR_UNLIKELY (MPFR_IS_NEG (a)))
     94     {
     95       MPFR_SET_NAN (r);
     96       MPFR_RET_NAN;
     97     }
     98 
     99   exp_a = MPFR_GET_EXP (a);
    100 
    101   /* If a is 1, the result is +0 */
    102   if (MPFR_UNLIKELY (exp_a == 1 && mpfr_cmp_ui (a, 1) == 0))
    103     {
    104       MPFR_SET_ZERO (r);
    105       MPFR_SET_POS (r);
    106       MPFR_RET (0); /* only "normal" case where the result is exact */
    107     }
    108 
    109   q = MPFR_PREC (r);
    110 
    111   /* use initial precision about q+2*lg(q)+cte */
    112   p = q + 2 * MPFR_INT_CEIL_LOG2 (q) + 10;
    113   /* % ~(mpfr_prec_t)GMP_NUMB_BITS  ;
    114      m=q; while (m) { p++; m >>= 1; }  */
    115   /* if (MPFR_LIKELY(p % GMP_NUMB_BITS != 0))
    116       p += GMP_NUMB_BITS - (p%GMP_NUMB_BITS); */
    117 
    118   MPFR_SAVE_EXPO_MARK (expo);
    119   MPFR_GROUP_INIT_2 (group, p, tmp1, tmp2);
    120 
    121   MPFR_ZIV_INIT (loop, p);
    122   for (;;)
    123     {
    124       mpfr_t scaled_a;
    125       mpfr_exp_t m;
    126       mpfr_exp_t cancel;
    127 
    128       /* Calculus of m (depends on p)
    129          If mpfr_exp_t has N bits, then both (p + 3) / 2 and |exp_a| fit
    130          on N-2 bits, so that there cannot be an overflow. */
    131       m = (p + 3) / 2 - exp_a;
    132 
    133       /* In standard configuration (_MPFR_EXP_FORMAT <= 3), one has
    134          mpfr_exp_t <= long, so that the following assertion is always
    135          true. This assertion is needed for the mpfr_mul_si below. */
    136       MPFR_ASSERTN (m >= LONG_MIN && m <= LONG_MAX);
    137 
    138       /* FIXME: Redo the error analysis. The error concerning the AGM
    139          should be explained since 4/s is inexact (one needs a bound
    140          on its derivative). */
    141       MPFR_ALIAS (scaled_a, a, MPFR_SIGN_POS, (p + 3) / 2); /* s=a*2^m */
    142       /* [FIXME] and one can have the equality, even if p is even.
    143          This means that if a is a power of 2 and p is even, then
    144          s = (1/2) * 2^((p+2)/2) = 2^(p/2), so that the condition
    145          s > 2^(p/2) from algorithms.tex is not satisfied. */
    146       mpfr_div (tmp1, __gmpfr_four, scaled_a, MPFR_RNDF); /* 4/s, err<=2 ulps */
    147       mpfr_agm (tmp2, __gmpfr_one, tmp1, MPFR_RNDN); /* AG(1,4/s),err<=3 ulps */
    148       mpfr_mul_2ui (tmp2, tmp2, 1, MPFR_RNDN); /* 2*AG(1,4/s),    err<=3 ulps */
    149       mpfr_const_pi (tmp1, MPFR_RNDN);         /* compute pi,     err<=1ulp   */
    150       mpfr_div (tmp2, tmp1, tmp2, MPFR_RNDN);  /* pi/2*AG(1,4/s), err<=5ulps  */
    151       mpfr_const_log2 (tmp1, MPFR_RNDN);      /* compute log(2),  err<=1ulp   */
    152       mpfr_mul_si (tmp1, tmp1, m, MPFR_RNDN); /* compute m*log(2),err<=2ulps  */
    153       mpfr_sub (tmp1, tmp2, tmp1, MPFR_RNDN); /* log(a),    err<=7ulps+cancel */
    154 
    155       if (MPFR_LIKELY (MPFR_IS_PURE_FP (tmp1) && MPFR_IS_PURE_FP (tmp2)))
    156         {
    157           cancel = MPFR_GET_EXP (tmp2) - MPFR_GET_EXP (tmp1);
    158           MPFR_LOG_MSG (("canceled bits=%ld\n", (long) cancel));
    159           MPFR_LOG_VAR (tmp1);
    160           if (MPFR_UNLIKELY (cancel < 0))
    161             cancel = 0;
    162 
    163           /* we have 7 ulps of error from the above roundings,
    164              4 ulps from the 4/s^2 second order term,
    165              plus the canceled bits */
    166           if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp1, p - cancel - 4, q, rnd_mode)))
    167             break;
    168 
    169           /* VL: I think it is better to have an increment that it isn't
    170              too low; in particular, the increment must be positive even
    171              if cancel = 0 (can this occur?). */
    172           p += cancel + MPFR_INT_CEIL_LOG2 (p);
    173         }
    174       else
    175         {
    176           /* TODO: find why this case can occur and what is best to do
    177              with it. */
    178           p += MPFR_INT_CEIL_LOG2 (p);
    179         }
    180 
    181       MPFR_ZIV_NEXT (loop, p);
    182       MPFR_GROUP_REPREC_2 (group, p, tmp1, tmp2);
    183     }
    184   MPFR_ZIV_FREE (loop);
    185   inexact = mpfr_set (r, tmp1, rnd_mode);
    186   /* We clean */
    187   MPFR_GROUP_CLEAR (group);
    188 
    189   MPFR_SAVE_EXPO_FREE (expo);
    190   return mpfr_check_range (r, inexact, rnd_mode);
    191 }
    192