set_uj.c revision 1.1.1.3
1/* mpfr_set_uj -- set a MPFR number from a huge machine unsigned integer 2 3Copyright 2004-2018 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 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" 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 33#define uintmaxpml (sizeof(uintmax_t) / sizeof(mp_limb_t)) 34 35int 36mpfr_set_uj (mpfr_t x, uintmax_t j, mpfr_rnd_t rnd) 37{ 38 return mpfr_set_uj_2exp (x, j, 0, rnd); 39} 40 41int 42mpfr_set_uj_2exp (mpfr_t x, uintmax_t j, intmax_t e, mpfr_rnd_t rnd) 43{ 44 int cnt; 45 mp_size_t i, k; 46 mp_limb_t limb; 47 mp_limb_t yp[uintmaxpml]; 48 mpfr_t y; 49 unsigned long uintmax_bit_size = sizeof(uintmax_t) * CHAR_BIT; 50 unsigned long bpml = GMP_NUMB_BITS % uintmax_bit_size; 51 52 /* Special case */ 53 if (j == 0) 54 { 55 MPFR_SET_POS(x); 56 MPFR_SET_ZERO(x); 57 MPFR_RET(0); 58 } 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 /* Check expo underflow / overflow (can't use mpfr_check_range) */ 111 if (MPFR_UNLIKELY(e < __gmpfr_emin)) 112 { 113 /* The following test is necessary because in the rounding to the 114 * nearest mode, mpfr_underflow always rounds away from 0. In 115 * this rounding mode, we need to round to 0 if: 116 * _ |x| < 2^(emin-2), or 117 * _ |x| = 2^(emin-2) and the absolute value of the exact 118 * result is <= 2^(emin-2). */ 119 if (rnd == MPFR_RNDN && (e+1 < __gmpfr_emin || mpfr_powerof2_raw(y))) 120 rnd = MPFR_RNDZ; 121 return mpfr_underflow (x, rnd, MPFR_SIGN_POS); 122 } 123 if (MPFR_UNLIKELY(e > __gmpfr_emax)) 124 return mpfr_overflow (x, rnd, MPFR_SIGN_POS); 125 MPFR_SET_EXP (y, e); 126 127 /* Final: set x to y (rounding if necessary) */ 128 return mpfr_set (x, y, rnd); 129} 130 131#endif 132