1/* mpfr_set_uj -- set a MPFR number from a huge machine unsigned integer
2
3Copyright 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
4Contributed by the Arenaire and Cacao 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/* The ISO C99 standard specifies that in C++ implementations the
28   INTMAX_MAX, ... macros should only be defined if explicitly requested.  */
29#if defined(__cplusplus)
30# define __STDC_LIMIT_MACROS
31# define __STDC_CONSTANT_MACROS
32#endif
33
34#if HAVE_INTTYPES_H
35# include <inttypes.h> /* for intmax_t */
36#else
37# if HAVE_STDINT_H
38#  include <stdint.h>
39# endif
40#endif
41
42#define MPFR_NEED_LONGLONG_H
43#include "mpfr-impl.h"
44
45#ifdef _MPFR_H_HAVE_INTMAX_T
46
47int
48mpfr_set_uj (mpfr_t x, uintmax_t j, mpfr_rnd_t rnd)
49{
50  return mpfr_set_uj_2exp (x, j, 0, rnd);
51}
52
53int
54mpfr_set_uj_2exp (mpfr_t x, uintmax_t j, intmax_t e, mpfr_rnd_t rnd)
55{
56  unsigned int cnt, i;
57  mp_size_t k, len;
58  mp_limb_t limb;
59  mp_limb_t yp[sizeof(uintmax_t) / sizeof(mp_limb_t)];
60  mpfr_t y;
61  unsigned long uintmax_bit_size = sizeof(uintmax_t) * CHAR_BIT;
62  unsigned long bpml = GMP_NUMB_BITS % uintmax_bit_size;
63
64  /* Special case */
65  if (j == 0)
66    {
67      MPFR_SET_POS(x);
68      MPFR_SET_ZERO(x);
69      MPFR_RET(0);
70    }
71
72  MPFR_ASSERTN (sizeof(uintmax_t) % sizeof(mp_limb_t) == 0);
73
74  /* Create an auxillary var */
75  MPFR_TMP_INIT1 (yp, y, uintmax_bit_size);
76  k = numberof (yp);
77  if (k == 1)
78    limb = yp[0] = j;
79  else
80    {
81      /* Note: either GMP_NUMB_BITS = uintmax_bit_size, then k = 1 the
82         shift j >>= bpml is never done, or GMP_NUMB_BITS < uintmax_bit_size
83         and bpml = GMP_NUMB_BITS. */
84      for (i = 0; i < k; i++, j >>= bpml)
85        yp[i] = j; /* Only the low bits are copied */
86
87      /* Find the first limb not equal to zero. */
88      do
89        {
90          MPFR_ASSERTD (k > 0);
91          limb = yp[--k];
92        }
93      while (limb == 0);
94      k++;
95    }
96  count_leading_zeros(cnt, limb);
97  len = numberof (yp) - k;
98
99  /* Normalize it: len = number of last 0 limb, k number of non-zero limbs */
100  if (MPFR_LIKELY(cnt))
101    mpn_lshift (yp+len, yp, k, cnt);  /* Normalize the High Limb*/
102  else if (len != 0)
103    MPN_COPY_DECR (yp+len, yp, k);    /* Must use DECR */
104  if (len != 0)
105    /* Note: when numberof(yp)==1, len is constant and null, so the compiler
106       can optimize out this code. */
107    {
108      if (len == 1)
109        yp[0] = (mp_limb_t) 0;
110      else
111        MPN_ZERO (yp, len);   /* Zeroing the last limbs */
112    }
113  e += k * GMP_NUMB_BITS - cnt;    /* Update Expo */
114  MPFR_ASSERTD (MPFR_LIMB_MSB(yp[numberof (yp) - 1]) != 0);
115
116  /* Check expo underflow / overflow (can't use mpfr_check_range) */
117  if (MPFR_UNLIKELY(e < __gmpfr_emin))
118    {
119      /* The following test is necessary because in the rounding to the
120       * nearest mode, mpfr_underflow always rounds away from 0. In
121       * this rounding mode, we need to round to 0 if:
122       *   _ |x| < 2^(emin-2), or
123       *   _ |x| = 2^(emin-2) and the absolute value of the exact
124       *     result is <= 2^(emin-2). */
125      if (rnd == MPFR_RNDN && (e+1 < __gmpfr_emin || mpfr_powerof2_raw(y)))
126        rnd = MPFR_RNDZ;
127      return mpfr_underflow (x, rnd, MPFR_SIGN_POS);
128    }
129  if (MPFR_UNLIKELY(e > __gmpfr_emax))
130    return mpfr_overflow (x, rnd, MPFR_SIGN_POS);
131  MPFR_SET_EXP (y, e);
132
133  /* Final: set x to y (rounding if necessary) */
134  return mpfr_set (x, y, rnd);
135}
136
137#endif
138