1 /* mpfr_mul_ui -- multiply a floating-point number by a machine integer
2
3 Copyright 1999-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 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
27 MPFR_HOT_FUNCTION_ATTR int
28 mpfr_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 }