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[%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 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 } 259