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