Index: gcc/mpfr/const_log2.c |
diff --git a/gcc/mpfr/const_log2.c b/gcc/mpfr/const_log2.c |
deleted file mode 100644 |
index 198ac6d677f0bd4c2b65b9c77c3425d14df0bd41..0000000000000000000000000000000000000000 |
--- a/gcc/mpfr/const_log2.c |
+++ /dev/null |
@@ -1,192 +0,0 @@ |
-/* mpfr_const_log2 -- compute natural logarithm of 2 |
- |
-Copyright 1999, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc. |
-Contributed by the Arenaire and Cacao projects, INRIA. |
- |
-This file is part of the GNU MPFR Library. |
- |
-The GNU MPFR Library is free software; you can redistribute it and/or modify |
-it under the terms of the GNU Lesser General Public License as published by |
-the Free Software Foundation; either version 2.1 of the License, or (at your |
-option) any later version. |
- |
-The GNU MPFR Library is distributed in the hope that it will be useful, but |
-WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
-or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public |
-License for more details. |
- |
-You should have received a copy of the GNU Lesser General Public License |
-along with the GNU MPFR Library; see the file COPYING.LIB. If not, write to |
-the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, |
-MA 02110-1301, USA. */ |
- |
-#define MPFR_NEED_LONGLONG_H |
-#include "mpfr-impl.h" |
- |
-/* Declare the cache */ |
-MPFR_DECL_INIT_CACHE(__gmpfr_cache_const_log2, mpfr_const_log2_internal); |
- |
-/* Set User interface */ |
-#undef mpfr_const_log2 |
-int |
-mpfr_const_log2 (mpfr_ptr x, mp_rnd_t rnd_mode) { |
- return mpfr_cache (x, __gmpfr_cache_const_log2, rnd_mode); |
-} |
- |
-/* Auxiliary function: Compute the terms from n1 to n2 (excluded) |
- 3/4*sum((-1)^n*n!^2/2^n/(2*n+1)!, n = n1..n2-1). |
- Numerator is T[0], denominator is Q[0], |
- Compute P[0] only when need_P is non-zero. |
- Need 1+ceil(log(n2-n1)/log(2)) cells in T[],P[],Q[]. |
-*/ |
-static void |
-S (mpz_t *T, mpz_t *P, mpz_t *Q, unsigned long n1, unsigned long n2, int need_P) |
-{ |
- if (n2 == n1 + 1) |
- { |
- if (n1 == 0) |
- mpz_set_ui (P[0], 3); |
- else |
- { |
- mpz_set_ui (P[0], n1); |
- mpz_neg (P[0], P[0]); |
- } |
- if (n1 <= (ULONG_MAX / 4 - 1) / 2) |
- mpz_set_ui (Q[0], 4 * (2 * n1 + 1)); |
- else /* to avoid overflow in 4 * (2 * n1 + 1) */ |
- { |
- mpz_set_ui (Q[0], n1); |
- mpz_mul_2exp (Q[0], Q[0], 1); |
- mpz_add_ui (Q[0], Q[0], 1); |
- mpz_mul_2exp (Q[0], Q[0], 2); |
- } |
- mpz_set (T[0], P[0]); |
- } |
- else |
- { |
- unsigned long m = (n1 / 2) + (n2 / 2) + (n1 & 1UL & n2); |
- unsigned long v, w; |
- |
- S (T, P, Q, n1, m, 1); |
- S (T + 1, P + 1, Q + 1, m, n2, need_P); |
- mpz_mul (T[0], T[0], Q[1]); |
- mpz_mul (T[1], T[1], P[0]); |
- mpz_add (T[0], T[0], T[1]); |
- if (need_P) |
- mpz_mul (P[0], P[0], P[1]); |
- mpz_mul (Q[0], Q[0], Q[1]); |
- |
- /* remove common trailing zeroes if any */ |
- v = mpz_scan1 (T[0], 0); |
- if (v > 0) |
- { |
- w = mpz_scan1 (Q[0], 0); |
- if (w < v) |
- v = w; |
- if (need_P) |
- { |
- w = mpz_scan1 (P[0], 0); |
- if (w < v) |
- v = w; |
- } |
- /* now v = min(val(T), val(Q), val(P)) */ |
- if (v > 0) |
- { |
- mpz_div_2exp (T[0], T[0], v); |
- mpz_div_2exp (Q[0], Q[0], v); |
- if (need_P) |
- mpz_div_2exp (P[0], P[0], v); |
- } |
- } |
- } |
-} |
- |
-/* Don't need to save / restore exponent range: the cache does it */ |
-int |
-mpfr_const_log2_internal (mpfr_ptr x, mp_rnd_t rnd_mode) |
-{ |
- unsigned long n = MPFR_PREC (x); |
- mp_prec_t w; /* working precision */ |
- unsigned long N; |
- mpz_t *T, *P, *Q; |
- mpfr_t t, q; |
- int inexact; |
- int ok = 1; /* ensures that the 1st try will give correct rounding */ |
- unsigned long lgN, i; |
- MPFR_ZIV_DECL (loop); |
- |
- MPFR_LOG_FUNC (("rnd_mode=%d", rnd_mode), ("x[%#R]=%R inex=%d",x,x,inexact)); |
- |
- mpfr_init2 (t, MPFR_PREC_MIN); |
- mpfr_init2 (q, MPFR_PREC_MIN); |
- |
- if (n < 1253) |
- w = n + 10; /* ensures correct rounding for the four rounding modes, |
- together with N = w / 3 + 1 (see below). */ |
- else if (n < 2571) |
- w = n + 11; /* idem */ |
- else if (n < 3983) |
- w = n + 12; |
- else if (n < 4854) |
- w = n + 13; |
- else if (n < 26248) |
- w = n + 14; |
- else |
- { |
- w = n + 15; |
- ok = 0; |
- } |
- |
- MPFR_ZIV_INIT (loop, w); |
- for (;;) |
- { |
- N = w / 3 + 1; /* Warning: do not change that (even increasing N!) |
- without checking correct rounding in the above |
- ranges for n. */ |
- |
- /* the following are needed for error analysis (see algorithms.tex) */ |
- MPFR_ASSERTD(w >= 3 && N >= 2); |
- |
- lgN = MPFR_INT_CEIL_LOG2 (N) + 1; |
- T = (mpz_t *) (*__gmp_allocate_func) (3 * lgN * sizeof (mpz_t)); |
- P = T + lgN; |
- Q = T + 2*lgN; |
- for (i = 0; i < lgN; i++) |
- { |
- mpz_init (T[i]); |
- mpz_init (P[i]); |
- mpz_init (Q[i]); |
- } |
- |
- S (T, P, Q, 0, N, 0); |
- |
- mpfr_set_prec (t, w); |
- mpfr_set_prec (q, w); |
- |
- mpfr_set_z (t, T[0], GMP_RNDN); |
- mpfr_set_z (q, Q[0], GMP_RNDN); |
- mpfr_div (t, t, q, GMP_RNDN); |
- |
- for (i = 0; i < lgN; i++) |
- { |
- mpz_clear (T[i]); |
- mpz_clear (P[i]); |
- mpz_clear (Q[i]); |
- } |
- (*__gmp_free_func) (T, 3 * lgN * sizeof (mpz_t)); |
- |
- if (MPFR_LIKELY (ok != 0 |
- || mpfr_can_round (t, w - 2, GMP_RNDN, rnd_mode, n))) |
- break; |
- |
- MPFR_ZIV_NEXT (loop, w); |
- } |
- MPFR_ZIV_FREE (loop); |
- |
- inexact = mpfr_set (x, t, rnd_mode); |
- |
- mpfr_clear (t); |
- mpfr_clear (q); |
- |
- return inexact; |
-} |