1 /* mpfr_const_log2 -- compute natural logarithm of 2
3 Copyright 1999, 2001-2015 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 /* Declare the cache */
27 #ifndef MPFR_USE_LOGGING
28 MPFR_DECL_INIT_CACHE(__gmpfr_cache_const_log2
, mpfr_const_log2_internal
);
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
;
35 /* Set User interface */
36 #undef mpfr_const_log2
38 mpfr_const_log2 (mpfr_ptr x
, mpfr_rnd_t rnd_mode
) {
39 return mpfr_cache (x
, __gmpfr_cache_const_log2
, rnd_mode
);
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[].
49 S (mpz_t
*T
, mpz_t
*P
, mpz_t
*Q
, unsigned long n1
, unsigned long n2
, int need_P
)
57 mpz_set_ui (P
[0], n1
);
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) */
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);
73 unsigned long m
= (n1
/ 2) + (n2
/ 2) + (n1
& 1UL & n2
);
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]);
82 mpz_mul (P
[0], P
[0], P
[1]);
83 mpz_mul (Q
[0], Q
[0], Q
[1]);
85 /* remove common trailing zeroes if any */
86 v
= mpz_scan1 (T
[0], 0);
89 w
= mpz_scan1 (Q
[0], 0);
94 w
= mpz_scan1 (P
[0], 0);
98 /* now v = min(val(T), val(Q), val(P)) */
101 mpz_fdiv_q_2exp (T
[0], T
[0], v
);
102 mpz_fdiv_q_2exp (Q
[0], Q
[0], v
);
104 mpz_fdiv_q_2exp (P
[0], P
[0], v
);
110 /* Don't need to save / restore exponent range: the cache does it */
112 mpfr_const_log2_internal (mpfr_ptr x
, mpfr_rnd_t rnd_mode
)
114 unsigned long n
= MPFR_PREC (x
);
115 mpfr_prec_t w
; /* working precision */
120 int ok
= 1; /* ensures that the 1st try will give correct rounding */
121 unsigned long lgN
, i
;
122 MPFR_ZIV_DECL (loop
);
125 ("rnd_mode=%d", rnd_mode
),
126 ("x[%Pu]=%.*Rg inex=%d", mpfr_get_prec(x
), mpfr_log_prec
, x
, inexact
));
128 mpfr_init2 (t
, MPFR_PREC_MIN
);
129 mpfr_init2 (q
, MPFR_PREC_MIN
);
132 w
= n
+ 10; /* ensures correct rounding for the four rounding modes,
133 together with N = w / 3 + 1 (see below). */
135 w
= n
+ 11; /* idem */
148 MPFR_ZIV_INIT (loop
, w
);
151 N
= w
/ 3 + 1; /* Warning: do not change that (even increasing N!)
152 without checking correct rounding in the above
155 /* the following are needed for error analysis (see algorithms.tex) */
156 MPFR_ASSERTD(w
>= 3 && N
>= 2);
158 lgN
= MPFR_INT_CEIL_LOG2 (N
) + 1;
159 T
= (mpz_t
*) (*__gmp_allocate_func
) (3 * lgN
* sizeof (mpz_t
));
162 for (i
= 0; i
< lgN
; i
++)
169 S (T
, P
, Q
, 0, N
, 0);
171 mpfr_set_prec (t
, w
);
172 mpfr_set_prec (q
, w
);
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
);
178 for (i
= 0; i
< lgN
; i
++)
184 (*__gmp_free_func
) (T
, 3 * lgN
* sizeof (mpz_t
));
186 if (MPFR_LIKELY (ok
!= 0
187 || mpfr_can_round (t
, w
- 2, MPFR_RNDN
, rnd_mode
, n
)))
190 MPFR_ZIV_NEXT (loop
, w
);
192 MPFR_ZIV_FREE (loop
);
194 inexact
= mpfr_set (x
, t
, rnd_mode
);