contrib/mpfr/src/get_f.c

   1 /* mpfr_get_f -- convert a MPFR number to a GNU MPF number
   2
   3 Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
   4 Contributed by the AriC and Caramel 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 3 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.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. */
  22
  23 #include "mpfr-impl.h"
  24
  25 /* Since MPFR-3.0, return the usual inexact value.
  26    The erange flag is set if an error occurred in the conversion
  27    (y is NaN, +Inf, or -Inf that have no equivalent in mpf)
  28 */
  29 int
  30 mpfr_get_f (mpf_ptr x, mpfr_srcptr y, mpfr_rnd_t rnd_mode)
  31 {
  32   int inex;
  33   mp_size_t sx, sy;
  34   mpfr_prec_t precx, precy;
  35   mp_limb_t *xp;
  36   int sh;
  37
  38   if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(y)))
  39     {
  40       if (MPFR_IS_ZERO(y))
  41         {
  42           mpf_set_ui (x, 0);
  43           return 0;
  44         }
  45       else if (MPFR_IS_NAN (y))
  46         {
  47           MPFR_SET_ERANGE ();
  48           return 0;
  49         }
  50       else /* y is plus infinity (resp. minus infinity), set x to the maximum
  51               value (resp. the minimum value) in precision PREC(x) */
  52         {
  53           int i;
  54           mp_limb_t *xp;
  55
  56           MPFR_SET_ERANGE ();
  57
  58           /* To this day, [mp_exp_t] and mp_size_t are #defined as the same
  59              type */
  60           EXP (x) = MP_SIZE_T_MAX;
  61
  62           sx = PREC (x);
  63           SIZ (x) = sx;
  64           xp = PTR (x);
  65           for (i = 0; i < sx; i++)
  66             xp[i] = MP_LIMB_T_MAX;
  67
  68           if (MPFR_IS_POS (y))
  69             return -1;
  70           else
  71             {
  72               mpf_neg (x, x);
  73               return +1;
  74             }
  75         }
  76     }
  77
  78   sx = PREC(x); /* number of limbs of the mantissa of x */
  79
  80   precy = MPFR_PREC(y);
  81   precx = (mpfr_prec_t) sx * GMP_NUMB_BITS;
  82   sy = MPFR_LIMB_SIZE (y);
  83
  84   xp = PTR (x);
  85
  86   /* since mpf numbers are represented in base 2^GMP_NUMB_BITS,
  87      we loose -EXP(y) % GMP_NUMB_BITS bits in the most significant limb */
  88   sh = MPFR_GET_EXP(y) % GMP_NUMB_BITS;
  89   sh = sh <= 0 ? - sh : GMP_NUMB_BITS - sh;
  90   MPFR_ASSERTD (sh >= 0);
  91   if (precy + sh <= precx) /* we can copy directly */
  92     {
  93       mp_size_t ds;
  94
  95       MPFR_ASSERTN (sx >= sy);
  96       ds = sx - sy;
  97
  98       if (sh != 0)
  99         {
 100           mp_limb_t out;
 101           out = mpn_rshift (xp + ds, MPFR_MANT(y), sy, sh);
 102           MPFR_ASSERTN (ds > 0 || out == 0);
 103           if (ds > 0)
 104             xp[--ds] = out;
 105         }
 106       else
 107         MPN_COPY (xp + ds, MPFR_MANT (y), sy);
 108       if (ds > 0)
 109         MPN_ZERO (xp, ds);
 110       EXP(x) = (MPFR_GET_EXP(y) + sh) / GMP_NUMB_BITS;
 111       inex = 0;
 112     }
 113   else /* we have to round to precx - sh bits */
 114     {
 115       mpfr_t z;
 116       mp_size_t sz;
 117
 118       /* Recall that precx = (mpfr_prec_t) sx * GMP_NUMB_BITS, thus removing
 119          sh bits (sh < GMP_NUMB_BITSS) won't reduce the number of limbs. */
 120       mpfr_init2 (z, precx - sh);
 121       sz = MPFR_LIMB_SIZE (z);
 122       MPFR_ASSERTN (sx == sz);
 123
 124       inex = mpfr_set (z, y, rnd_mode);
 125       /* warning, sh may change due to rounding, but then z is a power of two,
 126          thus we can safely ignore its last bit which is 0 */
 127       sh = MPFR_GET_EXP(z) % GMP_NUMB_BITS;
 128       sh = sh <= 0 ? - sh : GMP_NUMB_BITS - sh;
 129       MPFR_ASSERTD (sh >= 0);
 130       if (sh != 0)
 131         {
 132           mp_limb_t out;
 133           out = mpn_rshift (xp, MPFR_MANT(z), sz, sh);
 134           /* If sh hasn't changed, it is the number of the non-significant
 135              bits in the lowest limb of z. Therefore out == 0. */
 136           MPFR_ASSERTD (out == 0);  (void) out; /* avoid a warning */
 137         }
 138       else
 139         MPN_COPY (xp, MPFR_MANT(z), sz);
 140       EXP(x) = (MPFR_GET_EXP(z) + sh) / GMP_NUMB_BITS;
 141       mpfr_clear (z);
 142     }
 143
 144   /* set size and sign */
 145   SIZ(x) = (MPFR_FROM_SIGN_TO_INT(MPFR_SIGN(y)) < 0) ? -sx : sx;
 146
 147   return inex;
 148 }