1 /* Functions to work with unbounded floats (limited low-level interface).
2
3 Copyright 2016-2020 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 /* Note: In MPFR math functions, even if UBF code is not called first,
27 we may still need to handle special values NaN and infinities here.
28 Indeed, for MAXR * MAXR + (-inf), even though this is a special case,
29 the computation will also generate an overflow due to MAXR * MAXR,
30 so that UBF code will be called anyway (except if special cases are
31 detected and handled separately, but for polynomial, this should not
32 be needed). */
33
34 /* Get the exponent of a regular MPFR number or UBF as a mpz_t, which is
35 initialized by this function. The flags are not changed. */
36 static void
mpfr_init_get_zexp(mpz_ptr ez,mpfr_srcptr x)37 mpfr_init_get_zexp (mpz_ptr ez, mpfr_srcptr x)
38 {
39 mpz_init (ez);
40
41 if (MPFR_IS_UBF (x))
42 mpz_set (ez, MPFR_ZEXP (x));
43 else
44 {
45 mpfr_exp_t ex = MPFR_GET_EXP (x);
46
47 #if _MPFR_EXP_FORMAT <= 3
48 /* mpfr_exp_t fits in a long */
49 mpz_set_si (ez, ex);
50 #else
51 mp_limb_t e_limb[MPFR_EXP_LIMB_SIZE];
52 mpfr_t e;
53 int inex;
54 MPFR_SAVE_EXPO_DECL (expo);
55
56 MPFR_TMP_INIT1 (e_limb, e, sizeof (mpfr_exp_t) * CHAR_BIT);
57 MPFR_SAVE_EXPO_MARK (expo);
58 MPFR_DBGRES (inex = mpfr_set_exp_t (e, ex, MPFR_RNDN));
59 MPFR_ASSERTD (inex == 0);
60 MPFR_DBGRES (inex = mpfr_get_z (ez, e, MPFR_RNDN));
61 MPFR_ASSERTD (inex == 0);
62 MPFR_SAVE_EXPO_FREE (expo);
63 #endif
64 }
65 }
66
67 /* Exact product. The number a is assumed to have enough allocated memory,
68 where the trailing bits are regarded as being part of the input numbers
69 (no reallocation is attempted and no check is performed as MPFR_TMP_INIT
70 could have been used). The arguments b and c may actually be UBF numbers
71 (mpfr_srcptr can be seen a bit like void *, but is stronger).
72 This function does not change the flags, except in case of NaN. */
73 void
mpfr_ubf_mul_exact(mpfr_ubf_ptr a,mpfr_srcptr b,mpfr_srcptr c)74 mpfr_ubf_mul_exact (mpfr_ubf_ptr a, mpfr_srcptr b, mpfr_srcptr c)
75 {
76 MPFR_LOG_FUNC
77 (("b[%Pu]=%.*Rg c[%Pu]=%.*Rg",
78 mpfr_get_prec (b), mpfr_log_prec, b,
79 mpfr_get_prec (c), mpfr_log_prec, c),
80 ("a[%Pu]=%.*Rg",
81 mpfr_get_prec ((mpfr_ptr) a), mpfr_log_prec, a));
82
83 MPFR_ASSERTD ((mpfr_ptr) a != b);
84 MPFR_ASSERTD ((mpfr_ptr) a != c);
85 MPFR_SIGN (a) = MPFR_MULT_SIGN (MPFR_SIGN (b), MPFR_SIGN (c));
86
87 if (MPFR_ARE_SINGULAR (b, c))
88 {
89 if (MPFR_IS_NAN (b) || MPFR_IS_NAN (c))
90 MPFR_SET_NAN (a);
91 else if (MPFR_IS_INF (b))
92 {
93 if (MPFR_NOTZERO (c))
94 MPFR_SET_INF (a);
95 else
96 MPFR_SET_NAN (a);
97 }
98 else if (MPFR_IS_INF (c))
99 {
100 if (!MPFR_IS_ZERO (b))
101 MPFR_SET_INF (a);
102 else
103 MPFR_SET_NAN (a);
104 }
105 else
106 {
107 MPFR_ASSERTD (MPFR_IS_ZERO(b) || MPFR_IS_ZERO(c));
108 MPFR_SET_ZERO (a);
109 }
110 }
111 else
112 {
113 mpfr_exp_t e;
114 mp_size_t bn, cn;
115 mpfr_limb_ptr ap;
116 mp_limb_t u, v;
117 int m;
118
119 /* Note about the code below: For the choice of the precision of
120 * the result a, one could choose PREC(b) + PREC(c), instead of
121 * taking whole limbs into account, but in most cases where one
122 * would gain one limb, one would need to copy the significand
123 * instead of a no-op (see the mul.c code).
124 * But in the case MPFR_LIMB_MSB (u) == 0, if the result fits in
125 * an-1 limbs, one could actually do
126 * mpn_rshift (ap, ap, k, GMP_NUMB_BITS - 1)
127 * instead of
128 * mpn_lshift (ap, ap, k, 1)
129 * to gain one limb (and reduce the precision), replacing a shift
130 * by another one. Would this be interesting?
131 */
132
133 bn = MPFR_LIMB_SIZE (b);
134 cn = MPFR_LIMB_SIZE (c);
135
136 ap = MPFR_MANT (a);
137
138 if (bn == 1 && cn == 1)
139 {
140 umul_ppmm (ap[1], ap[0], MPFR_MANT(b)[0], MPFR_MANT(c)[0]);
141 if (ap[1] & MPFR_LIMB_HIGHBIT)
142 m = 0;
143 else
144 {
145 ap[1] = (ap[1] << 1) | (ap[0] >> (GMP_NUMB_BITS - 1));
146 ap[0] = ap[0] << 1;
147 m = 1;
148 }
149 }
150 else
151 {
152 u = (bn >= cn) ?
153 mpn_mul (ap, MPFR_MANT (b), bn, MPFR_MANT (c), cn) :
154 mpn_mul (ap, MPFR_MANT (c), cn, MPFR_MANT (b), bn);
155 if (MPFR_LIMB_MSB (u) == 0)
156 {
157 m = 1;
158 MPFR_DBGRES (v = mpn_lshift (ap, ap, bn + cn, 1));
159 MPFR_ASSERTD (v == 0);
160 }
161 else
162 m = 0;
163 }
164
165 if (! MPFR_IS_UBF (b) && ! MPFR_IS_UBF (c) &&
166 (e = MPFR_GET_EXP (b) + MPFR_GET_EXP (c) - m,
167 MPFR_EXP_IN_RANGE (e)))
168 {
169 MPFR_SET_EXP (a, e);
170 }
171 else
172 {
173 mpz_t be, ce;
174
175 mpz_init (MPFR_ZEXP (a));
176
177 /* This may involve copies of mpz_t, but exponents should not be
178 very large integers anyway. */
179 mpfr_init_get_zexp (be, b);
180 mpfr_init_get_zexp (ce, c);
181 mpz_add (MPFR_ZEXP (a), be, ce);
182 mpz_clear (be);
183 mpz_clear (ce);
184 mpz_sub_ui (MPFR_ZEXP (a), MPFR_ZEXP (a), m);
185 MPFR_SET_UBF (a);
186 }
187 }
188 }
189
190 /* Compare the exponents of two numbers, which can be either MPFR numbers
191 or UBF numbers. If both numbers can be MPFR numbers, it is better to
192 use the MPFR_UBF_EXP_LESS_P wrapper macro, which is optimized for this
193 common case. */
194 int
mpfr_ubf_exp_less_p(mpfr_srcptr x,mpfr_srcptr y)195 mpfr_ubf_exp_less_p (mpfr_srcptr x, mpfr_srcptr y)
196 {
197 mpz_t xe, ye;
198 int c;
199
200 mpfr_init_get_zexp (xe, x);
201 mpfr_init_get_zexp (ye, y);
202 c = mpz_cmp (xe, ye) < 0;
203 mpz_clear (xe);
204 mpz_clear (ye);
205 return c;
206 }
207
208 /* Convert an mpz_t to an mpfr_exp_t, saturated to
209 the interval [MPFR_EXP_MIN,MPFR_EXP_MAX]. */
210 mpfr_exp_t
mpfr_ubf_zexp2exp(mpz_ptr ez)211 mpfr_ubf_zexp2exp (mpz_ptr ez)
212 {
213 mp_size_t n;
214 mpfr_eexp_t e;
215 mpfr_t d;
216 int inex;
217 MPFR_SAVE_EXPO_DECL (expo);
218
219 n = ABSIZ (ez); /* limb size of ez */
220 if (n == 0)
221 return 0;
222
223 MPFR_SAVE_EXPO_MARK (expo);
224 mpfr_init2 (d, n * GMP_NUMB_BITS);
225 MPFR_DBGRES (inex = mpfr_set_z (d, ez, MPFR_RNDN));
226 MPFR_ASSERTD (inex == 0);
227 e = mpfr_get_exp_t (d, MPFR_RNDZ);
228 mpfr_clear (d);
229 MPFR_SAVE_EXPO_FREE (expo);
230 if (MPFR_UNLIKELY (e < MPFR_EXP_MIN))
231 return MPFR_EXP_MIN;
232 if (MPFR_UNLIKELY (e > MPFR_EXP_MAX))
233 return MPFR_EXP_MAX;
234 return e;
235 }
236
237 /* Return the difference of the exponents of x and y, saturated to
238 the interval [MPFR_EXP_MIN,MPFR_EXP_MAX]. */
239 mpfr_exp_t
mpfr_ubf_diff_exp(mpfr_srcptr x,mpfr_srcptr y)240 mpfr_ubf_diff_exp (mpfr_srcptr x, mpfr_srcptr y)
241 {
242 mpz_t xe, ye;
243 mpfr_exp_t e;
244
245 mpfr_init_get_zexp (xe, x);
246 mpfr_init_get_zexp (ye, y);
247 mpz_sub (xe, xe, ye);
248 mpz_clear (ye);
249 e = mpfr_ubf_zexp2exp (xe);
250 mpz_clear (xe);
251 return e;
252 }
253