1 /* Functions to work with unbounded floats (limited low-level interface).
2
3 Copyright 2016-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 /* 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
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
74 mpfr_ubf_mul_exact (mpfr_ubf_ptr a, mpfr_srcptr b, mpfr_srcptr c)
75 {
76 MPFR_LOG_FUNC
77 (("b[%Pd]=%.*Rg c[%Pd]=%.*Rg",
78 mpfr_get_prec (b), mpfr_log_prec, b,
79 mpfr_get_prec (c), mpfr_log_prec, c),
80 ("a[%Pd]=%.*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 if (b == c)
153 {
154 mpn_sqr (ap, MPFR_MANT (b), bn);
155 u = ap[2 * bn - 1];
156 }
157 else
158 u = (bn >= cn) ?
159 mpn_mul (ap, MPFR_MANT (b), bn, MPFR_MANT (c), cn) :
160 mpn_mul (ap, MPFR_MANT (c), cn, MPFR_MANT (b), bn);
161 if (MPFR_LIMB_MSB (u) == 0)
162 {
163 m = 1;
164 MPFR_DBGRES (v = mpn_lshift (ap, ap, bn + cn, 1));
165 MPFR_ASSERTD (v == 0);
166 }
167 else
168 m = 0;
169 }
170
171 if (! MPFR_IS_UBF (b) && ! MPFR_IS_UBF (c) &&
172 (e = MPFR_GET_EXP (b) + MPFR_GET_EXP (c) - m,
173 MPFR_EXP_IN_RANGE (e)))
174 {
175 MPFR_SET_EXP (a, e);
176 }
177 else
178 {
179 mpz_t be, ce;
180
181 mpz_init (MPFR_ZEXP (a));
182
183 /* This may involve copies of mpz_t, but exponents should not be
184 very large integers anyway. */
185 mpfr_init_get_zexp (be, b);
186 mpfr_init_get_zexp (ce, c);
187 mpz_add (MPFR_ZEXP (a), be, ce);
188 mpz_clear (be);
189 mpz_clear (ce);
190 mpz_sub_ui (MPFR_ZEXP (a), MPFR_ZEXP (a), m);
191 MPFR_SET_UBF (a);
192 }
193 }
194 }
195
196 /* Compare the exponents of two numbers, which can be either MPFR numbers
197 or UBF numbers. If both numbers can be MPFR numbers, it is better to
198 use the MPFR_UBF_EXP_LESS_P wrapper macro, which is optimized for this
199 common case. */
200 int
201 mpfr_ubf_exp_less_p (mpfr_srcptr x, mpfr_srcptr y)
202 {
203 mpz_t xe, ye;
204 int c;
205
206 mpfr_init_get_zexp (xe, x);
207 mpfr_init_get_zexp (ye, y);
208 c = mpz_cmp (xe, ye) < 0;
209 mpz_clear (xe);
210 mpz_clear (ye);
211 return c;
212 }
213
214 /* Convert an mpz_t to an mpfr_exp_t, saturated to
215 the interval [MPFR_EXP_MIN,MPFR_EXP_MAX]. */
216 mpfr_exp_t
217 mpfr_ubf_zexp2exp (mpz_ptr ez)
218 {
219 mp_size_t n;
220 mpfr_eexp_t e;
221 mpfr_t d;
222 int inex;
223 MPFR_SAVE_EXPO_DECL (expo);
224
225 n = ABSIZ (ez); /* limb size of ez */
226 if (n == 0)
227 return 0;
228
229 MPFR_SAVE_EXPO_MARK (expo);
230 mpfr_init2 (d, n * GMP_NUMB_BITS);
231 MPFR_DBGRES (inex = mpfr_set_z (d, ez, MPFR_RNDN));
232 MPFR_ASSERTD (inex == 0);
233 e = mpfr_get_exp_t (d, MPFR_RNDZ);
234 mpfr_clear (d);
235 MPFR_SAVE_EXPO_FREE (expo);
236 if (MPFR_UNLIKELY (e < MPFR_EXP_MIN))
237 return MPFR_EXP_MIN;
238 if (MPFR_UNLIKELY (e > MPFR_EXP_MAX))
239 return MPFR_EXP_MAX;
240 return e;
241 }
242
243 /* Return the difference of the exponents of x and y, saturated to
244 the interval [MPFR_EXP_MIN,MPFR_EXP_MAX]. */
245 mpfr_exp_t
246 mpfr_ubf_diff_exp (mpfr_srcptr x, mpfr_srcptr y)
247 {
248 mpz_t xe, ye;
249 mpfr_exp_t e;
250
251 mpfr_init_get_zexp (xe, x);
252 mpfr_init_get_zexp (ye, y);
253 mpz_sub (xe, xe, ye);
254 mpz_clear (ye);
255 e = mpfr_ubf_zexp2exp (xe);
256 mpz_clear (xe);
257 return e;
258 }