mul_ui.c revision 1.1.1.4
1/* mpfr_mul_ui -- multiply a floating-point number by a machine integer 2 3Copyright 1999-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#include "mpfr-impl.h" 25 26#undef mpfr_mul_ui 27MPFR_HOT_FUNCTION_ATTR int 28mpfr_mul_ui (mpfr_ptr y, mpfr_srcptr x, unsigned long int u, mpfr_rnd_t rnd_mode) 29{ 30 int inexact; 31 32 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 33 { 34 if (MPFR_IS_NAN (x)) 35 { 36 MPFR_SET_NAN (y); 37 MPFR_RET_NAN; 38 } 39 else if (MPFR_IS_INF (x)) 40 { 41 if (u != 0) 42 { 43 MPFR_SET_INF (y); 44 MPFR_SET_SAME_SIGN (y, x); 45 MPFR_RET (0); /* infinity is exact */ 46 } 47 else /* 0 * infinity */ 48 { 49 MPFR_SET_NAN (y); 50 MPFR_RET_NAN; 51 } 52 } 53 else /* x is zero */ 54 { 55 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 56 MPFR_SET_ZERO (y); 57 MPFR_SET_SAME_SIGN (y, x); 58 MPFR_RET (0); /* zero is exact */ 59 } 60 } 61 else if (u <= 1) 62 { 63 if (u < 1) 64 { 65 MPFR_SET_ZERO (y); 66 MPFR_SET_SAME_SIGN (y, x); 67 MPFR_RET (0); /* zero is exact */ 68 } 69 else 70 return mpfr_set (y, x, rnd_mode); 71 } 72 else if (MPFR_UNLIKELY (IS_POW2 (u))) 73 return mpfr_mul_2si (y, x, MPFR_INT_CEIL_LOG2 (u), rnd_mode); 74 75#ifdef MPFR_LONG_WITHIN_LIMB 76 { 77 mp_limb_t *yp; 78 mp_size_t xn; 79 int cnt; 80 MPFR_TMP_DECL (marker); 81 82 yp = MPFR_MANT (y); 83 xn = MPFR_LIMB_SIZE (x); 84 85 MPFR_ASSERTD (xn < MP_SIZE_T_MAX); 86 MPFR_TMP_MARK(marker); 87 yp = MPFR_TMP_LIMBS_ALLOC (xn + 1); 88 89 MPFR_ASSERTN (u == (mp_limb_t) u); 90 yp[xn] = mpn_mul_1 (yp, MPFR_MANT (x), xn, u); 91 92 /* x * u is stored in yp[xn], ..., yp[0] */ 93 94 /* since the case u=1 was treated above, we have u >= 2, thus 95 yp[xn] >= 1 since x was msb-normalized */ 96 MPFR_ASSERTD (yp[xn] != 0); 97 if (MPFR_LIKELY (MPFR_LIMB_MSB (yp[xn]) == 0)) 98 { 99 count_leading_zeros (cnt, yp[xn]); 100 mpn_lshift (yp, yp, xn + 1, cnt); 101 } 102 else 103 { 104 cnt = 0; 105 } 106 107 /* now yp[xn], ..., yp[0] is msb-normalized too, and has at most 108 PREC(x) + (GMP_NUMB_BITS - cnt) non-zero bits */ 109 MPFR_RNDRAW (inexact, y, yp, (mpfr_prec_t) (xn + 1) * GMP_NUMB_BITS, 110 rnd_mode, MPFR_SIGN (x), cnt -- ); 111 112 MPFR_TMP_FREE (marker); 113 114 cnt = GMP_NUMB_BITS - cnt; 115 if (MPFR_UNLIKELY (__gmpfr_emax < MPFR_EMAX_MIN + cnt 116 || MPFR_GET_EXP (x) > __gmpfr_emax - cnt)) 117 return mpfr_overflow (y, rnd_mode, MPFR_SIGN(x)); 118 119 MPFR_SET_EXP (y, MPFR_GET_EXP (x) + cnt); 120 MPFR_SET_SAME_SIGN (y, x); 121 MPFR_RET (inexact); 122 } 123#else 124 { 125 mpfr_t uu; 126 MPFR_SAVE_EXPO_DECL (expo); 127 128 mpfr_init2 (uu, sizeof (unsigned long) * CHAR_BIT); 129 /* Warning: u might be outside the current exponent range! */ 130 MPFR_SAVE_EXPO_MARK (expo); 131 mpfr_set_ui (uu, u, MPFR_RNDZ); 132 inexact = mpfr_mul (y, x, uu, rnd_mode); 133 mpfr_clear (uu); 134 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); 135 MPFR_SAVE_EXPO_FREE (expo); 136 return mpfr_check_range (y, inexact, rnd_mode); 137 } 138#endif /* MPFR_LONG_WITHIN_LIMB */ 139} 140