contrib/mpfr/set_z.c

   1 /* mpfr_set_z -- set a floating-point number from a multiple-precision integer
   2
   3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao 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 2.1 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.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 MA 02110-1301, USA. */
  22
  23 #define MPFR_NEED_LONGLONG_H
  24 #include "mpfr-impl.h"
  25
  26 /* set f to the integer z */
  27 int
  28 mpfr_set_z (mpfr_ptr f, mpz_srcptr z, mp_rnd_t rnd_mode)
  29 {
  30   mp_size_t fn, zn, dif;
  31   int k, sign_z, inex;
  32   mp_limb_t *fp, *zp;
  33   mp_exp_t exp;
  34
  35   sign_z = mpz_sgn (z);
  36   if (MPFR_UNLIKELY (sign_z == 0))
  37     {
  38       MPFR_SET_ZERO(f);
  39       MPFR_SET_POS(f);
  40       MPFR_RET(0);
  41     }
  42   MPFR_ASSERTD (sign_z == MPFR_SIGN_POS || sign_z == MPFR_SIGN_NEG);
  43
  44   zn = ABS(SIZ(z));
  45   MPFR_ASSERTD (zn >= 1);
  46   if (MPFR_UNLIKELY (zn > MPFR_EMAX_MAX / BITS_PER_MP_LIMB + 1))
  47     return mpfr_overflow (f, rnd_mode, sign_z);
  48
  49   fp = MPFR_MANT (f);
  50   fn = MPFR_LIMB_SIZE (f);
  51   dif = zn - fn;
  52   zp = PTR(z);
  53   count_leading_zeros (k, zp[zn-1]);
  54   /* because zn >= __gmpfr_emax / BITS_PER_MP_LIMB + 2
  55      and zn * BITS_PER_MP_LIMB >= __gmpfr_emax + BITS_PER_MP_LIMB + 1
  56      and exp = zn * BITS_PER_MP_LIMB - k > __gmpfr_emax */
  57
  58   /* now zn <= MPFR_EMAX_MAX / BITS_PER_MP_LIMB + 1
  59      thus zn * BITS_PER_MP_LIMB <= MPFR_EMAX_MAX + BITS_PER_MP_LIMB
  60      and exp = zn * BITS_PER_MP_LIMB - k
  61              <= MPFR_EMAX_MAX + BITS_PER_MP_LIMB */
  62   exp = (mp_prec_t) zn * BITS_PER_MP_LIMB - k;
  63   /* The exponent will be exp or exp + 1 (due to rounding) */
  64   if (MPFR_UNLIKELY (exp > __gmpfr_emax))
  65     return mpfr_overflow (f, rnd_mode, sign_z);
  66   if (MPFR_UNLIKELY (exp + 1 < __gmpfr_emin))
  67     return mpfr_underflow(f, rnd_mode == GMP_RNDN ? GMP_RNDZ : rnd_mode,
  68                               sign_z);
  69
  70   if (MPFR_LIKELY (dif >= 0))
  71     {
  72       mp_limb_t rb, sb, ulp;
  73       int sh;
  74
  75       /* number has to be truncated */
  76       if (MPFR_LIKELY (k != 0))
  77         {
  78           mpn_lshift (fp, &zp[dif], fn, k);
  79           if (MPFR_LIKELY (dif > 0))
  80             fp[0] |= zp[dif - 1] >> (BITS_PER_MP_LIMB - k);
  81         }
  82       else
  83         MPN_COPY (fp, zp + dif, fn);
  84
  85       /* Compute Rounding Bit and Sticky Bit */
  86       MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (f) );
  87       if (MPFR_LIKELY (sh != 0))
  88         {
  89           mp_limb_t mask = MPFR_LIMB_ONE << (sh-1);
  90           mp_limb_t limb = fp[0];
  91           rb = limb & mask;
  92           sb = limb & (mask-1);
  93           ulp = 2*mask;
  94           fp[0] = limb & ~(ulp-1);
  95         }
  96       else /* sh == 0 */
  97         {
  98           mp_limb_t mask = MPFR_LIMB_ONE << (BITS_PER_MP_LIMB - 1 - k);
  99           if (MPFR_LIKELY (dif > 0))
 100             {
 101               rb = zp[--dif] & mask;
 102               sb = zp[dif] & (mask-1);
 103             }
 104           else
 105             rb = sb = 0;
 106           k = 0;
 107           ulp = MPFR_LIMB_ONE;
 108         }
 109       if (MPFR_UNLIKELY (sb == 0) && MPFR_LIKELY (dif > 0))
 110         {
 111           sb = zp[--dif];
 112           if (MPFR_LIKELY (k != 0))
 113             sb &= MPFR_LIMB_MASK (BITS_PER_MP_LIMB - k);
 114           if (MPFR_UNLIKELY (sb == 0) && MPFR_LIKELY (dif > 0))
 115             do {
 116               sb = zp[--dif];
 117             } while (dif > 0 && sb == 0);
 118         }
 119
 120       /* Rounding */
 121       if (MPFR_LIKELY (rnd_mode == GMP_RNDN))
 122         {
 123           if (rb == 0 || MPFR_UNLIKELY (sb == 0 && (fp[0] & ulp) == 0))
 124             goto trunc;
 125           else
 126             goto addoneulp;
 127         }
 128       else /* Not Nearest */
 129         {
 130           if (MPFR_LIKELY (MPFR_IS_LIKE_RNDZ (rnd_mode, sign_z < 0))
 131               || MPFR_UNLIKELY ( (sb|rb) == 0 ))
 132             goto trunc;
 133           else
 134             goto addoneulp;
 135         }
 136
 137     trunc:
 138       inex = MPFR_LIKELY ((sb | rb) != 0) ? -1 : 0;
 139       goto end;
 140
 141     addoneulp:
 142       inex = 1;
 143       if (MPFR_UNLIKELY (mpn_add_1 (fp, fp, fn, ulp)))
 144         {
 145           /* Pow 2 case */
 146           if (MPFR_UNLIKELY (exp == __gmpfr_emax))
 147             return mpfr_overflow (f, rnd_mode, sign_z);
 148           exp ++;
 149           fp[fn-1] = MPFR_LIMB_HIGHBIT;
 150         }
 151     end:
 152       (void) 0;
 153     }
 154   else   /* dif < 0: Mantissa F is strictly bigger than z's one */
 155     {
 156       if (MPFR_LIKELY (k != 0))
 157         mpn_lshift (fp - dif, zp, zn, k);
 158       else
 159         MPN_COPY (fp - dif, zp, zn);
 160       /* fill with zeroes */
 161       MPN_ZERO (fp, -dif);
 162       inex = 0; /* result is exact */
 163     }
 164
 165   if (MPFR_UNLIKELY (exp < __gmpfr_emin))
 166     {
 167       if (rnd_mode == GMP_RNDN && inex == 0 && mpfr_powerof2_raw (f))
 168         rnd_mode = GMP_RNDZ;
 169       return mpfr_underflow(f, rnd_mode, sign_z);
 170     }
 171
 172   MPFR_SET_EXP (f, exp);
 173   MPFR_SET_SIGN (f, sign_z);
 174   MPFR_RET (inex*sign_z);
 175 }
 176