1/* mpfr_const_log2 -- compute natural logarithm of 2
2
3Copyright 1999, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
4Contributed by the Arenaire and Cacao projects, INRIA.
5
6This file is part of the GNU MPFR Library.
7
8The GNU MPFR Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MPFR Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2151 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 */
27MPFR_DECL_INIT_CACHE(__gmpfr_cache_const_log2, mpfr_const_log2_internal);
28
29/* Set User interface */
30#undef mpfr_const_log2
31int
32mpfr_const_log2 (mpfr_ptr x, mpfr_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*/
42static void
43S (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_fdiv_q_2exp (T[0], T[0], v);
96              mpz_fdiv_q_2exp (Q[0], Q[0], v);
97              if (need_P)
98                mpz_fdiv_q_2exp (P[0], P[0], v);
99            }
100        }
101    }
102}
103
104/* Don't need to save / restore exponent range: the cache does it */
105int
106mpfr_const_log2_internal (mpfr_ptr x, mpfr_rnd_t rnd_mode)
107{
108  unsigned long n = MPFR_PREC (x);
109  mpfr_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], MPFR_RNDN);
167      mpfr_set_z (q, Q[0], MPFR_RNDN);
168      mpfr_div (t, t, q, MPFR_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, MPFR_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}
193