| Index: gcc/mpfr/cbrt.c
|
| diff --git a/gcc/mpfr/cbrt.c b/gcc/mpfr/cbrt.c
|
| deleted file mode 100644
|
| index 923f4c304da64b57c9edd777bc56f4e97bdd60d1..0000000000000000000000000000000000000000
|
| --- a/gcc/mpfr/cbrt.c
|
| +++ /dev/null
|
| @@ -1,148 +0,0 @@
|
| -/* mpfr_cbrt -- cube root function.
|
| -
|
| -Copyright 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"
|
| -
|
| - /* The computation of y = x^(1/3) is done as follows:
|
| -
|
| - Let x = sign * m * 2^(3*e) where m is an integer
|
| -
|
| - with 2^(3n-3) <= m < 2^(3n) where n = PREC(y)
|
| -
|
| - and m = s^3 + r where 0 <= r and m < (s+1)^3
|
| -
|
| - we want that s has n bits i.e. s >= 2^(n-1), or m >= 2^(3n-3)
|
| - i.e. m must have at least 3n-2 bits
|
| -
|
| - then x^(1/3) = s * 2^e if r=0
|
| - x^(1/3) = (s+1) * 2^e if round up
|
| - x^(1/3) = (s-1) * 2^e if round down
|
| - x^(1/3) = s * 2^e if nearest and r < 3/2*s^2+3/4*s+1/8
|
| - (s+1) * 2^e otherwise
|
| - */
|
| -
|
| -int
|
| -mpfr_cbrt (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
|
| -{
|
| - mpz_t m;
|
| - mp_exp_t e, r, sh;
|
| - mp_prec_t n, size_m, tmp;
|
| - int inexact, negative;
|
| - MPFR_SAVE_EXPO_DECL (expo);
|
| -
|
| - /* special values */
|
| - if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
|
| - {
|
| - if (MPFR_IS_NAN (x))
|
| - {
|
| - MPFR_SET_NAN (y);
|
| - MPFR_RET_NAN;
|
| - }
|
| - else if (MPFR_IS_INF (x))
|
| - {
|
| - MPFR_SET_INF (y);
|
| - MPFR_SET_SAME_SIGN (y, x);
|
| - MPFR_RET (0);
|
| - }
|
| - /* case 0: cbrt(+/- 0) = +/- 0 */
|
| - else /* x is necessarily 0 */
|
| - {
|
| - MPFR_ASSERTD (MPFR_IS_ZERO (x));
|
| - MPFR_SET_ZERO (y);
|
| - MPFR_SET_SAME_SIGN (y, x);
|
| - MPFR_RET (0);
|
| - }
|
| - }
|
| -
|
| - /* General case */
|
| - MPFR_SAVE_EXPO_MARK (expo);
|
| - mpz_init (m);
|
| -
|
| - e = mpfr_get_z_exp (m, x); /* x = m * 2^e */
|
| - if ((negative = MPFR_IS_NEG(x)))
|
| - mpz_neg (m, m);
|
| - r = e % 3;
|
| - if (r < 0)
|
| - r += 3;
|
| - /* x = (m*2^r) * 2^(e-r) = (m*2^r) * 2^(3*q) */
|
| -
|
| - MPFR_MPZ_SIZEINBASE2 (size_m, m);
|
| - n = MPFR_PREC (y) + (rnd_mode == GMP_RNDN);
|
| -
|
| - /* we want 3*n-2 <= size_m + 3*sh + r <= 3*n
|
| - i.e. 3*sh + size_m + r <= 3*n */
|
| - sh = (3 * (mp_exp_t) n - (mp_exp_t) size_m - r) / 3;
|
| - sh = 3 * sh + r;
|
| - if (sh >= 0)
|
| - {
|
| - mpz_mul_2exp (m, m, sh);
|
| - e = e - sh;
|
| - }
|
| - else if (r > 0)
|
| - {
|
| - mpz_mul_2exp (m, m, r);
|
| - e = e - r;
|
| - }
|
| -
|
| - /* invariant: x = m*2^e, with e divisible by 3 */
|
| -
|
| - /* we reuse the variable m to store the cube root, since it is not needed
|
| - any more: we just need to know if the root is exact */
|
| - inexact = mpz_root (m, m, 3) == 0;
|
| -
|
| - MPFR_MPZ_SIZEINBASE2 (tmp, m);
|
| - sh = tmp - n;
|
| - if (sh > 0) /* we have to flush to 0 the last sh bits from m */
|
| - {
|
| - inexact = inexact || ((mp_exp_t) mpz_scan1 (m, 0) < sh);
|
| - mpz_div_2exp (m, m, sh);
|
| - e += 3 * sh;
|
| - }
|
| -
|
| - if (inexact)
|
| - {
|
| - if (negative)
|
| - rnd_mode = MPFR_INVERT_RND (rnd_mode);
|
| - if (rnd_mode == GMP_RNDU
|
| - || (rnd_mode == GMP_RNDN && mpz_tstbit (m, 0)))
|
| - inexact = 1, mpz_add_ui (m, m, 1);
|
| - else
|
| - inexact = -1;
|
| - }
|
| -
|
| - /* either inexact is not zero, and the conversion is exact, i.e. inexact
|
| - is not changed; or inexact=0, and inexact is set only when
|
| - rnd_mode=GMP_RNDN and bit (n+1) from m is 1 */
|
| - inexact += mpfr_set_z (y, m, GMP_RNDN);
|
| - MPFR_SET_EXP (y, MPFR_GET_EXP (y) + e / 3);
|
| -
|
| - if (negative)
|
| - {
|
| - MPFR_CHANGE_SIGN (y);
|
| - inexact = -inexact;
|
| - }
|
| -
|
| - mpz_clear (m);
|
| - MPFR_SAVE_EXPO_FREE (expo);
|
| - return mpfr_check_range (y, inexact, rnd_mode);
|
| -}
|
|
|