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