1/* mpfr_set_uj -- set a MPFR number from a huge machine unsigned integer
2
3Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4Contributed by the AriC and Caramel projects, INRIA.
5
6This file is part of the GNU MPFR Library.
7
8The GNU MPFR Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MPFR Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23#ifdef HAVE_CONFIG_H
24# include "config.h"       /* for a build within gmp */
25#endif
26
27#define MPFR_NEED_LONGLONG_H
28#include "mpfr-intmax.h"
29#include "mpfr-impl.h"
30
31#ifdef _MPFR_H_HAVE_INTMAX_T
32
33int
34mpfr_set_uj (mpfr_t x, uintmax_t j, mpfr_rnd_t rnd)
35{
36  return mpfr_set_uj_2exp (x, j, 0, rnd);
37}
38
39int
40mpfr_set_uj_2exp (mpfr_t x, uintmax_t j, intmax_t e, mpfr_rnd_t rnd)
41{
42  unsigned int cnt, i;
43  mp_size_t k, len;
44  mp_limb_t limb;
45  mp_limb_t yp[sizeof(uintmax_t) / sizeof(mp_limb_t)];
46  mpfr_t y;
47  unsigned long uintmax_bit_size = sizeof(uintmax_t) * CHAR_BIT;
48  unsigned long bpml = GMP_NUMB_BITS % uintmax_bit_size;
49
50  /* Special case */
51  if (j == 0)
52    {
53      MPFR_SET_POS(x);
54      MPFR_SET_ZERO(x);
55      MPFR_RET(0);
56    }
57
58  MPFR_ASSERTN (sizeof(uintmax_t) % sizeof(mp_limb_t) == 0);
59
60  /* Create an auxillary var */
61  MPFR_TMP_INIT1 (yp, y, uintmax_bit_size);
62  k = numberof (yp);
63  if (k == 1)
64    limb = yp[0] = j;
65  else
66    {
67      /* Note: either GMP_NUMB_BITS = uintmax_bit_size, then k = 1 the
68         shift j >>= bpml is never done, or GMP_NUMB_BITS < uintmax_bit_size
69         and bpml = GMP_NUMB_BITS. */
70      for (i = 0; i < k; i++, j >>= bpml)
71        yp[i] = j; /* Only the low bits are copied */
72
73      /* Find the first limb not equal to zero. */
74      do
75        {
76          MPFR_ASSERTD (k > 0);
77          limb = yp[--k];
78        }
79      while (limb == 0);
80      k++;
81    }
82  count_leading_zeros(cnt, limb);
83  len = numberof (yp) - k;
84
85  /* Normalize it: len = number of last 0 limb, k number of non-zero limbs */
86  if (MPFR_LIKELY(cnt))
87    mpn_lshift (yp+len, yp, k, cnt);  /* Normalize the High Limb*/
88  else if (len != 0)
89    MPN_COPY_DECR (yp+len, yp, k);    /* Must use DECR */
90  if (len != 0)
91    /* Note: when numberof(yp)==1, len is constant and null, so the compiler
92       can optimize out this code. */
93    {
94      if (len == 1)
95        yp[0] = (mp_limb_t) 0;
96      else
97        MPN_ZERO (yp, len);   /* Zeroing the last limbs */
98    }
99  e += k * GMP_NUMB_BITS - cnt;    /* Update Expo */
100  MPFR_ASSERTD (MPFR_LIMB_MSB(yp[numberof (yp) - 1]) != 0);
101
102  /* Check expo underflow / overflow (can't use mpfr_check_range) */
103  if (MPFR_UNLIKELY(e < __gmpfr_emin))
104    {
105      /* The following test is necessary because in the rounding to the
106       * nearest mode, mpfr_underflow always rounds away from 0. In
107       * this rounding mode, we need to round to 0 if:
108       *   _ |x| < 2^(emin-2), or
109       *   _ |x| = 2^(emin-2) and the absolute value of the exact
110       *     result is <= 2^(emin-2). */
111      if (rnd == MPFR_RNDN && (e+1 < __gmpfr_emin || mpfr_powerof2_raw(y)))
112        rnd = MPFR_RNDZ;
113      return mpfr_underflow (x, rnd, MPFR_SIGN_POS);
114    }
115  if (MPFR_UNLIKELY(e > __gmpfr_emax))
116    return mpfr_overflow (x, rnd, MPFR_SIGN_POS);
117  MPFR_SET_EXP (y, e);
118
119  /* Final: set x to y (rounding if necessary) */
120  return mpfr_set (x, y, rnd);
121}
122
123#endif
124