contrib/mpfr/src/const_log2.c

   1 /* mpfr_const_log2 -- compute natural logarithm of 2
   2
   3 Copyright 1999, 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 /* Declare the cache */
  27 #ifndef MPFR_USE_LOGGING
  28 MPFR_DECL_INIT_CACHE(__gmpfr_cache_const_log2, mpfr_const_log2_internal);
  29 #else
  30 MPFR_DECL_INIT_CACHE(__gmpfr_normal_log2, mpfr_const_log2_internal);
  31 MPFR_DECL_INIT_CACHE(__gmpfr_logging_log2, mpfr_const_log2_internal);
  32 mpfr_cache_ptr MPFR_THREAD_ATTR __gmpfr_cache_const_log2 = __gmpfr_normal_log2;
  33 #endif
  34
  35 /* Set User interface */
  36 #undef mpfr_const_log2
  37 int
  38 mpfr_const_log2 (mpfr_ptr x, mpfr_rnd_t rnd_mode) {
  39   return mpfr_cache (x, __gmpfr_cache_const_log2, rnd_mode);
  40 }
  41
  42 /* Auxiliary function: Compute the terms from n1 to n2 (excluded)
  43    3/4*sum((-1)^n*n!^2/2^n/(2*n+1)!, n = n1..n2-1).
  44    Numerator is T[0], denominator is Q[0],
  45    Compute P[0] only when need_P is non-zero.
  46    Need 1+ceil(log(n2-n1)/log(2)) cells in T[],P[],Q[].
  47 */
  48 static void
  49 S (mpz_t *T, mpz_t *P, mpz_t *Q, unsigned long n1, unsigned long n2, int need_P)
  50 {
  51   if (n2 == n1 + 1)
  52     {
  53       if (n1 == 0)
  54         mpz_set_ui (P[0], 3);
  55       else
  56         {
  57           mpz_set_ui (P[0], n1);
  58           mpz_neg (P[0], P[0]);
  59         }
  60       if (n1 <= (ULONG_MAX / 4 - 1) / 2)
  61         mpz_set_ui (Q[0], 4 * (2 * n1 + 1));
  62       else /* to avoid overflow in 4 * (2 * n1 + 1) */
  63         {
  64           mpz_set_ui (Q[0], n1);
  65           mpz_mul_2exp (Q[0], Q[0], 1);
  66           mpz_add_ui (Q[0], Q[0], 1);
  67           mpz_mul_2exp (Q[0], Q[0], 2);
  68         }
  69       mpz_set (T[0], P[0]);
  70     }
  71   else
  72     {
  73       unsigned long m = (n1 / 2) + (n2 / 2) + (n1 & 1UL & n2);
  74       unsigned long v, w;
  75
  76       S (T, P, Q, n1, m, 1);
  77       S (T + 1, P + 1, Q + 1, m, n2, need_P);
  78       mpz_mul (T[0], T[0], Q[1]);
  79       mpz_mul (T[1], T[1], P[0]);
  80       mpz_add (T[0], T[0], T[1]);
  81       if (need_P)
  82         mpz_mul (P[0], P[0], P[1]);
  83       mpz_mul (Q[0], Q[0], Q[1]);
  84
  85       /* remove common trailing zeroes if any */
  86       v = mpz_scan1 (T[0], 0);
  87       if (v > 0)
  88         {
  89           w = mpz_scan1 (Q[0], 0);
  90           if (w < v)
  91             v = w;
  92           if (need_P)
  93             {
  94               w = mpz_scan1 (P[0], 0);
  95               if (w < v)
  96                 v = w;
  97             }
  98           /* now v = min(val(T), val(Q), val(P)) */
  99           if (v > 0)
 100             {
 101               mpz_fdiv_q_2exp (T[0], T[0], v);
 102               mpz_fdiv_q_2exp (Q[0], Q[0], v);
 103               if (need_P)
 104                 mpz_fdiv_q_2exp (P[0], P[0], v);
 105             }
 106         }
 107     }
 108 }
 109
 110 /* Don't need to save / restore exponent range: the cache does it */
 111 int
 112 mpfr_const_log2_internal (mpfr_ptr x, mpfr_rnd_t rnd_mode)
 113 {
 114   unsigned long n = MPFR_PREC (x);
 115   mpfr_prec_t w; /* working precision */
 116   unsigned long N;
 117   mpz_t *T, *P, *Q;
 118   mpfr_t t, q;
 119   int inexact;
 120   int ok = 1; /* ensures that the 1st try will give correct rounding */
 121   unsigned long lgN, i;
 122   MPFR_ZIV_DECL (loop);
 123
 124   MPFR_LOG_FUNC (
 125     ("rnd_mode=%d", rnd_mode),
 126     ("x[%Pu]=%.*Rg inex=%d", mpfr_get_prec(x), mpfr_log_prec, x, inexact));
 127
 128   mpfr_init2 (t, MPFR_PREC_MIN);
 129   mpfr_init2 (q, MPFR_PREC_MIN);
 130
 131   if (n < 1253)
 132     w = n + 10; /* ensures correct rounding for the four rounding modes,
 133                    together with N = w / 3 + 1 (see below). */
 134   else if (n < 2571)
 135     w = n + 11; /* idem */
 136   else if (n < 3983)
 137     w = n + 12;
 138   else if (n < 4854)
 139     w = n + 13;
 140   else if (n < 26248)
 141     w = n + 14;
 142   else
 143     {
 144       w = n + 15;
 145       ok = 0;
 146     }
 147
 148   MPFR_ZIV_INIT (loop, w);
 149   for (;;)
 150     {
 151       N = w / 3 + 1; /* Warning: do not change that (even increasing N!)
 152                         without checking correct rounding in the above
 153                         ranges for n. */
 154
 155       /* the following are needed for error analysis (see algorithms.tex) */
 156       MPFR_ASSERTD(w >= 3 && N >= 2);
 157
 158       lgN = MPFR_INT_CEIL_LOG2 (N) + 1;
 159       T  = (mpz_t *) (*__gmp_allocate_func) (3 * lgN * sizeof (mpz_t));
 160       P  = T + lgN;
 161       Q  = T + 2*lgN;
 162       for (i = 0; i < lgN; i++)
 163         {
 164           mpz_init (T[i]);
 165           mpz_init (P[i]);
 166           mpz_init (Q[i]);
 167         }
 168
 169       S (T, P, Q, 0, N, 0);
 170
 171       mpfr_set_prec (t, w);
 172       mpfr_set_prec (q, w);
 173
 174       mpfr_set_z (t, T[0], MPFR_RNDN);
 175       mpfr_set_z (q, Q[0], MPFR_RNDN);
 176       mpfr_div (t, t, q, MPFR_RNDN);
 177
 178       for (i = 0; i < lgN; i++)
 179         {
 180           mpz_clear (T[i]);
 181           mpz_clear (P[i]);
 182           mpz_clear (Q[i]);
 183         }
 184       (*__gmp_free_func) (T, 3 * lgN * sizeof (mpz_t));
 185
 186       if (MPFR_LIKELY (ok != 0
 187                        || mpfr_can_round (t, w - 2, MPFR_RNDN, rnd_mode, n)))
 188         break;
 189
 190       MPFR_ZIV_NEXT (loop, w);
 191     }
 192   MPFR_ZIV_FREE (loop);
 193
 194   inexact = mpfr_set (x, t, rnd_mode);
 195
 196   mpfr_clear (t);
 197   mpfr_clear (q);
 198
 199   return inexact;
 200 }