1/* mpfr_set_uj -- set a MPFR number from a huge machine unsigned integer
2
3Copyright 2004-2023 Free Software Foundation, Inc.
4Contributed by the AriC and Caramba 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
20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23#define MPFR_NEED_LONGLONG_H
24#define MPFR_NEED_INTMAX_H
25#include "mpfr-impl.h"
26
27#ifdef _MPFR_H_HAVE_INTMAX_T
28
29#define uintmaxpml (sizeof(uintmax_t) / sizeof(mp_limb_t))
30
31int
32mpfr_set_uj (mpfr_ptr x, uintmax_t j, mpfr_rnd_t rnd)
33{
34  return mpfr_set_uj_2exp (x, j, 0, rnd);
35}
36
37int
38mpfr_set_uj_2exp (mpfr_ptr x, uintmax_t j, intmax_t e, mpfr_rnd_t rnd)
39{
40  int cnt, inex;
41  mp_size_t i, k;
42  mp_limb_t limb;
43  mp_limb_t yp[uintmaxpml];
44  mpfr_t y;
45  unsigned long uintmax_bit_size = sizeof(uintmax_t) * CHAR_BIT;
46  unsigned long bpml = GMP_NUMB_BITS % uintmax_bit_size;
47
48  /* Special case */
49  if (j == 0)
50    {
51      MPFR_SET_POS(x);
52      MPFR_SET_ZERO(x);
53      MPFR_RET(0);
54    }
55
56  /* early overflow detection to avoid a possible integer overflow below */
57  if (MPFR_UNLIKELY(e >= __gmpfr_emax))
58    return mpfr_overflow (x, rnd, MPFR_SIGN_POS);
59
60  MPFR_ASSERTN (sizeof(uintmax_t) % sizeof(mp_limb_t) == 0);
61
62  /* Create an auxiliary var */
63  MPFR_TMP_INIT1 (yp, y, uintmax_bit_size);
64  /* The compiler will optimize the code by removing the useless branch. */
65  k = uintmaxpml;
66  if (uintmaxpml == 1)
67    {
68      limb = j;
69      count_leading_zeros(cnt, limb);
70      /* Normalize the most significant limb */
71      yp[0] = limb << cnt;
72    }
73  else
74    {
75      mp_size_t len;
76      /* Note: either GMP_NUMB_BITS = uintmax_bit_size, then k = 1 the
77         shift j >>= bpml is never done, or GMP_NUMB_BITS < uintmax_bit_size
78         and bpml = GMP_NUMB_BITS. */
79      for (i = 0; i < k; i++, j >>= bpml)
80        yp[i] = j; /* Only the low bits are copied */
81
82      /* Find the first limb not equal to zero. */
83      do
84        {
85          MPFR_ASSERTD (k > 0);
86          limb = yp[--k];
87        }
88      while (limb == 0);
89      k++;
90      len = numberof (yp) - k;
91      count_leading_zeros(cnt, limb);
92
93      /* Normalize it: len = number of last zero limbs,
94         k = number of previous limbs */
95      if (MPFR_LIKELY (cnt != 0))
96        mpn_lshift (yp+len, yp, k, cnt);  /* Normalize the high limb */
97      else if (len != 0)
98        mpn_copyd (yp+len, yp, k);    /* Must use copyd */
99      if (len != 0)
100        {
101          if (len == 1)
102            yp[0] = MPFR_LIMB_ZERO;
103          else
104            MPN_ZERO (yp, len);   /* Zero the last limbs */
105        }
106    }
107  e += k * GMP_NUMB_BITS - cnt;    /* Update Expo */
108  MPFR_ASSERTD (MPFR_LIMB_MSB(yp[numberof (yp) - 1]) != 0);
109
110  MPFR_RNDRAW (inex, x, yp, uintmax_bit_size, rnd, MPFR_SIGN_POS, e++);
111
112  /* Check expo underflow / overflow */
113  if (MPFR_UNLIKELY(e < __gmpfr_emin))
114    {
115      /* The following test is necessary because in the rounding to the
116       * nearest mode, mpfr_underflow always rounds away from 0. In
117       * this rounding mode, we need to round to 0 if:
118       *   _ |x| < 2^(emin-2), or
119       *   _ |x| = 2^(emin-2) and the absolute value of the exact
120       *     result is <= 2^(emin-2). */
121      if (rnd == MPFR_RNDN &&
122          (e + 1 < __gmpfr_emin ||
123           (mpfr_powerof2_raw (x) && inex >= 0)))
124        rnd = MPFR_RNDZ;
125      return mpfr_underflow (x, rnd, MPFR_SIGN_POS);
126    }
127  if (MPFR_UNLIKELY(e > __gmpfr_emax))
128    return mpfr_overflow (x, rnd, MPFR_SIGN_POS);
129
130  MPFR_SET_SIGN (x, MPFR_SIGN_POS);
131  MPFR_SET_EXP (x, e);
132  MPFR_RET (inex);
133}
134
135#endif
136