set_uj.c revision 1.1.1.4
1/* mpfr_set_uj -- set a MPFR number from a huge machine unsigned integer 2 3Copyright 2004-2020 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_t 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_t 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