1/* Implementations of operations between mpfr and mpz/mpq data 2 3Copyright 2001, 2003-2023 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/* TODO: for functions with mpz_srcptr, check whether mpz_fits_slong_p 27 is really useful in all cases. For instance, concerning the addition, 28 one now has mpz_t -> long -> unsigned long -> mpfr_t then mpfr_add 29 instead of mpz_t -> mpfr_t then mpfr_add. */ 30 31/* Init and set a mpfr_t with enough precision to store a mpz. 32 This function should be called in the extended exponent range. */ 33static void 34init_set_z (mpfr_ptr t, mpz_srcptr z) 35{ 36 mpfr_prec_t p; 37 int i; 38 39 if (mpz_size (z) <= 1) 40 p = GMP_NUMB_BITS; 41 else 42 MPFR_MPZ_SIZEINBASE2 (p, z); 43 mpfr_init2 (t, p); 44 i = mpfr_set_z (t, z, MPFR_RNDN); 45 /* Possible assertion failure in case of overflow. Such cases, 46 which imply that z is huge (if the function is called in 47 the extended exponent range), are currently not supported, 48 just like precisions around MPFR_PREC_MAX. */ 49 MPFR_ASSERTN (i == 0); (void) i; /* use i to avoid a warning */ 50} 51 52/* Init, set a mpfr_t with enough precision to store a mpz_t without round, 53 call the function, and clear the allocated mpfr_t */ 54static int 55foo (mpfr_ptr x, mpfr_srcptr y, mpz_srcptr z, mpfr_rnd_t r, 56 int (*f)(mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t)) 57{ 58 mpfr_t t; 59 int i; 60 MPFR_SAVE_EXPO_DECL (expo); 61 62 MPFR_SAVE_EXPO_MARK (expo); 63 init_set_z (t, z); /* There should be no exceptions. */ 64 i = (*f) (x, y, t, r); 65 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); 66 mpfr_clear (t); 67 MPFR_SAVE_EXPO_FREE (expo); 68 return mpfr_check_range (x, i, r); 69} 70 71static int 72foo2 (mpfr_ptr x, mpz_srcptr y, mpfr_srcptr z, mpfr_rnd_t r, 73 int (*f)(mpfr_ptr, mpfr_srcptr, mpfr_srcptr, mpfr_rnd_t)) 74{ 75 mpfr_t t; 76 int i; 77 MPFR_SAVE_EXPO_DECL (expo); 78 79 MPFR_SAVE_EXPO_MARK (expo); 80 init_set_z (t, y); /* There should be no exceptions. */ 81 i = (*f) (x, t, z, r); 82 MPFR_SAVE_EXPO_UPDATE_FLAGS (expo, __gmpfr_flags); 83 mpfr_clear (t); 84 MPFR_SAVE_EXPO_FREE (expo); 85 return mpfr_check_range (x, i, r); 86} 87 88int 89mpfr_mul_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t r) 90{ 91 if (mpz_fits_slong_p (z)) 92 return mpfr_mul_si (y, x, mpz_get_si (z), r); 93 else 94 return foo (y, x, z, r, mpfr_mul); 95} 96 97int 98mpfr_div_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t r) 99{ 100 if (mpz_fits_slong_p (z)) 101 return mpfr_div_si (y, x, mpz_get_si (z), r); 102 else 103 return foo (y, x, z, r, mpfr_div); 104} 105 106int 107mpfr_add_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t r) 108{ 109 if (mpz_fits_slong_p (z)) 110 return mpfr_add_si (y, x, mpz_get_si (z), r); 111 else 112 return foo (y, x, z, r, mpfr_add); 113} 114 115int 116mpfr_sub_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mpfr_rnd_t r) 117{ 118 if (mpz_fits_slong_p (z)) 119 return mpfr_sub_si (y, x, mpz_get_si (z), r); 120 else 121 return foo (y, x, z, r, mpfr_sub); 122} 123 124int 125mpfr_z_sub (mpfr_ptr y, mpz_srcptr x, mpfr_srcptr z, mpfr_rnd_t r) 126{ 127 if (mpz_fits_slong_p (x)) 128 return mpfr_si_sub (y, mpz_get_si (x), z, r); 129 else 130 return foo2 (y, x, z, r, mpfr_sub); 131} 132 133int 134mpfr_cmp_z (mpfr_srcptr x, mpz_srcptr z) 135{ 136 mpfr_t t; 137 int res; 138 mpfr_prec_t p; 139 mpfr_flags_t flags; 140 141 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 142 return mpfr_cmp_si (x, mpz_sgn (z)); 143 144 if (mpz_fits_slong_p (z)) 145 return mpfr_cmp_si (x, mpz_get_si (z)); 146 147 if (mpz_size (z) <= 1) 148 p = GMP_NUMB_BITS; 149 else 150 MPFR_MPZ_SIZEINBASE2 (p, z); 151 mpfr_init2 (t, p); 152 flags = __gmpfr_flags; 153 if (mpfr_set_z (t, z, MPFR_RNDN)) 154 { 155 /* overflow (t is an infinity) or underflow: z does not fit in the 156 current exponent range. 157 If overflow, then z is larger than the largest *integer* < +Inf 158 (if z > 0), thus we get t = +Inf (or -Inf), and the value of 159 mpfr_cmp (x, t) below is correct. 160 If underflow, then z is smaller than the smallest number > 0, 161 which is necessarily an integer, say xmin. 162 If z > xmin/2, then t is xmin, and we divide t by 2 to ensure t 163 is zero, and then the value of mpfr_cmp (x, t) below is correct. */ 164 mpfr_div_2ui (t, t, 2, MPFR_RNDZ); /* if underflow, set t to zero */ 165 __gmpfr_flags = flags; /* restore the flags */ 166 /* The real value of t (= z), which falls outside the exponent range, 167 has been replaced by an equivalent value for the comparison: zero 168 or an infinity. */ 169 } 170 res = mpfr_cmp (x, t); 171 mpfr_clear (t); 172 return res; 173} 174 175#ifndef MPFR_USE_MINI_GMP 176/* Compute y = RND(x*n/d), where n and d are mpz integers. 177 An integer 0 is assumed to have a positive sign. 178 This function is used by mpfr_mul_q and mpfr_div_q. 179 Note: the status of the rational 0/(-1) is not clear (if there is 180 a signed infinity, there should be a signed zero). But infinities 181 are not currently supported/documented in GMP, and if the rational 182 is canonicalized as it should be, the case 0/(-1) cannot occur. */ 183static int 184mpfr_muldiv_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr n, mpz_srcptr d, 185 mpfr_rnd_t rnd_mode) 186{ 187 if (MPFR_UNLIKELY (mpz_sgn (n) == 0)) 188 { 189 if (MPFR_UNLIKELY (mpz_sgn (d) == 0)) 190 MPFR_SET_NAN (y); 191 else 192 { 193 mpfr_mul_ui (y, x, 0, MPFR_RNDN); /* exact: +0, -0 or NaN */ 194 if (MPFR_UNLIKELY (mpz_sgn (d) < 0)) 195 MPFR_CHANGE_SIGN (y); 196 } 197 return 0; 198 } 199 else if (MPFR_UNLIKELY (mpz_sgn (d) == 0)) 200 { 201 mpfr_div_ui (y, x, 0, MPFR_RNDN); /* exact: +Inf, -Inf or NaN */ 202 if (MPFR_UNLIKELY (mpz_sgn (n) < 0)) 203 MPFR_CHANGE_SIGN (y); 204 return 0; 205 } 206 else 207 { 208 mpfr_prec_t p; 209 mpfr_t tmp; 210 int inexact; 211 MPFR_SAVE_EXPO_DECL (expo); 212 213 MPFR_SAVE_EXPO_MARK (expo); 214 215 /* With the current MPFR code, using mpfr_mul_z and mpfr_div_z 216 for the general case should be faster than doing everything 217 in mpn, mpz and/or mpq. MPFR_SAVE_EXPO_MARK could be avoided 218 here, but it would be more difficult to handle corner cases. */ 219 MPFR_MPZ_SIZEINBASE2 (p, n); 220 mpfr_init2 (tmp, MPFR_PREC (x) + p); 221 inexact = mpfr_mul_z (tmp, x, n, MPFR_RNDN); 222 /* Since |n| >= 1, an underflow is not possible. And the precision of 223 tmp has been chosen so that inexact != 0 iff there's an overflow. */ 224 if (MPFR_UNLIKELY (inexact != 0)) 225 { 226 mpfr_t x0; 227 mpfr_exp_t ex; 228 MPFR_BLOCK_DECL (flags); 229 230 /* intermediate overflow case */ 231 MPFR_ASSERTD (mpfr_inf_p (tmp)); 232 ex = MPFR_GET_EXP (x); /* x is a pure FP number */ 233 MPFR_ALIAS (x0, x, MPFR_SIGN(x), 0); /* x0 = x / 2^ex */ 234 MPFR_BLOCK (flags, 235 inexact = mpfr_mul_z (tmp, x0, n, MPFR_RNDN); 236 MPFR_ASSERTD (inexact == 0); 237 inexact = mpfr_div_z (y, tmp, d, rnd_mode); 238 /* Just in case the division underflows 239 (highly unlikely, not supported)... */ 240 MPFR_ASSERTN (!MPFR_BLOCK_EXCEP)); 241 MPFR_EXP (y) += ex; 242 /* Detect highly unlikely, not supported corner cases... */ 243 MPFR_ASSERTN (MPFR_EXP (y) >= __gmpfr_emin); 244 MPFR_ASSERTN (! MPFR_IS_SINGULAR (y)); 245 /* The potential overflow will be detected by mpfr_check_range. */ 246 } 247 else 248 inexact = mpfr_div_z (y, tmp, d, rnd_mode); 249 250 mpfr_clear (tmp); 251 252 MPFR_SAVE_EXPO_FREE (expo); 253 return mpfr_check_range (y, inexact, rnd_mode); 254 } 255} 256 257int 258mpfr_mul_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mpfr_rnd_t rnd_mode) 259{ 260 return mpfr_muldiv_z (y, x, mpq_numref (z), mpq_denref (z), rnd_mode); 261} 262 263int 264mpfr_div_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mpfr_rnd_t rnd_mode) 265{ 266 return mpfr_muldiv_z (y, x, mpq_denref (z), mpq_numref (z), rnd_mode); 267} 268 269int 270mpfr_add_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mpfr_rnd_t rnd_mode) 271{ 272 mpfr_t t,q; 273 mpfr_prec_t p; 274 mpfr_exp_t err; 275 int res; 276 MPFR_SAVE_EXPO_DECL (expo); 277 MPFR_ZIV_DECL (loop); 278 279 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 280 { 281 if (MPFR_IS_NAN (x)) 282 { 283 MPFR_SET_NAN (y); 284 MPFR_RET_NAN; 285 } 286 else if (MPFR_IS_INF (x)) 287 { 288 if (MPFR_UNLIKELY (mpz_sgn (mpq_denref (z)) == 0 && 289 MPFR_MULT_SIGN (mpz_sgn (mpq_numref (z)), 290 MPFR_SIGN (x)) <= 0)) 291 { 292 MPFR_SET_NAN (y); 293 MPFR_RET_NAN; 294 } 295 MPFR_SET_INF (y); 296 MPFR_SET_SAME_SIGN (y, x); 297 MPFR_RET (0); 298 } 299 else 300 { 301 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 302 if (MPFR_UNLIKELY (mpq_sgn (z) == 0)) 303 return mpfr_set (y, x, rnd_mode); /* signed 0 - Unsigned 0 */ 304 else 305 return mpfr_set_q (y, z, rnd_mode); 306 } 307 } 308 309 MPFR_SAVE_EXPO_MARK (expo); 310 311 p = MPFR_PREC (y) + 10; 312 mpfr_init2 (t, p); 313 mpfr_init2 (q, p); 314 315 MPFR_ZIV_INIT (loop, p); 316 for (;;) 317 { 318 MPFR_BLOCK_DECL (flags); 319 320 res = mpfr_set_q (q, z, MPFR_RNDN); /* Error <= 1/2 ulp(q) */ 321 /* If z if @INF@ (1/0), res = 0, so it quits immediately */ 322 if (MPFR_UNLIKELY (res == 0)) 323 /* Result is exact so we can add it directly! */ 324 { 325 res = mpfr_add (y, x, q, rnd_mode); 326 break; 327 } 328 MPFR_BLOCK (flags, mpfr_add (t, x, q, MPFR_RNDN)); 329 /* Error on t is <= 1/2 ulp(t), except in case of overflow/underflow, 330 but such an exception is very unlikely as it would be possible 331 only if q has a huge numerator or denominator. Not supported! */ 332 MPFR_ASSERTN (! (MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags))); 333 /* Error / ulp(t) <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t)) 334 If EXP(q)-EXP(t)>0, <= 2^(EXP(q)-EXP(t)-1)*(1+2^-(EXP(q)-EXP(t))) 335 <= 2^(EXP(q)-EXP(t)) 336 If EXP(q)-EXP(t)<0, <= 2^0 */ 337 /* We can get 0, but we can't round since q is inexact */ 338 if (MPFR_LIKELY (!MPFR_IS_ZERO (t))) 339 { 340 err = (mpfr_exp_t) p - 1 - MAX (MPFR_GET_EXP(q)-MPFR_GET_EXP(t), 0); 341 if (MPFR_LIKELY (MPFR_CAN_ROUND (t, err, MPFR_PREC (y), rnd_mode))) 342 { 343 res = mpfr_set (y, t, rnd_mode); 344 break; 345 } 346 } 347 MPFR_ZIV_NEXT (loop, p); 348 mpfr_set_prec (t, p); 349 mpfr_set_prec (q, p); 350 } 351 MPFR_ZIV_FREE (loop); 352 mpfr_clear (t); 353 mpfr_clear (q); 354 355 MPFR_SAVE_EXPO_FREE (expo); 356 return mpfr_check_range (y, res, rnd_mode); 357} 358 359int 360mpfr_sub_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z,mpfr_rnd_t rnd_mode) 361{ 362 mpfr_t t,q; 363 mpfr_prec_t p; 364 int res; 365 mpfr_exp_t err; 366 MPFR_SAVE_EXPO_DECL (expo); 367 MPFR_ZIV_DECL (loop); 368 369 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 370 { 371 if (MPFR_IS_NAN (x)) 372 { 373 MPFR_SET_NAN (y); 374 MPFR_RET_NAN; 375 } 376 else if (MPFR_IS_INF (x)) 377 { 378 if (MPFR_UNLIKELY (mpz_sgn (mpq_denref (z)) == 0 && 379 MPFR_MULT_SIGN (mpz_sgn (mpq_numref (z)), 380 MPFR_SIGN (x)) >= 0)) 381 { 382 MPFR_SET_NAN (y); 383 MPFR_RET_NAN; 384 } 385 MPFR_SET_INF (y); 386 MPFR_SET_SAME_SIGN (y, x); 387 MPFR_RET (0); 388 } 389 else 390 { 391 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 392 393 if (MPFR_UNLIKELY (mpq_sgn (z) == 0)) 394 return mpfr_set (y, x, rnd_mode); /* signed 0 - Unsigned 0 */ 395 else 396 { 397 res = mpfr_set_q (y, z, MPFR_INVERT_RND (rnd_mode)); 398 MPFR_CHANGE_SIGN (y); 399 return -res; 400 } 401 } 402 } 403 404 MPFR_SAVE_EXPO_MARK (expo); 405 406 p = MPFR_PREC (y) + 10; 407 mpfr_init2 (t, p); 408 mpfr_init2 (q, p); 409 410 MPFR_ZIV_INIT (loop, p); 411 for(;;) 412 { 413 MPFR_BLOCK_DECL (flags); 414 415 res = mpfr_set_q(q, z, MPFR_RNDN); /* Error <= 1/2 ulp(q) */ 416 /* If z if @INF@ (1/0), res = 0, so it quits immediately */ 417 if (MPFR_UNLIKELY (res == 0)) 418 /* Result is exact so we can add it directly!*/ 419 { 420 res = mpfr_sub (y, x, q, rnd_mode); 421 break; 422 } 423 MPFR_BLOCK (flags, mpfr_sub (t, x, q, MPFR_RNDN)); 424 /* Error on t is <= 1/2 ulp(t), except in case of overflow/underflow, 425 but such an exception is very unlikely as it would be possible 426 only if q has a huge numerator or denominator. Not supported! */ 427 MPFR_ASSERTN (! (MPFR_OVERFLOW (flags) || MPFR_UNDERFLOW (flags))); 428 /* Error / ulp(t) <= 1/2 + 1/2 * 2^(EXP(q)-EXP(t)) 429 If EXP(q)-EXP(t)>0, <= 2^(EXP(q)-EXP(t)-1)*(1+2^-(EXP(q)-EXP(t))) 430 <= 2^(EXP(q)-EXP(t)) 431 If EXP(q)-EXP(t)<0, <= 2^0 */ 432 /* We can get 0, but we can't round since q is inexact */ 433 if (MPFR_LIKELY (!MPFR_IS_ZERO (t))) 434 { 435 err = (mpfr_exp_t) p - 1 - MAX (MPFR_GET_EXP(q)-MPFR_GET_EXP(t), 0); 436 res = MPFR_CAN_ROUND (t, err, MPFR_PREC (y), rnd_mode); 437 if (MPFR_LIKELY (res != 0)) /* We can round! */ 438 { 439 res = mpfr_set (y, t, rnd_mode); 440 break; 441 } 442 } 443 MPFR_ZIV_NEXT (loop, p); 444 mpfr_set_prec (t, p); 445 mpfr_set_prec (q, p); 446 } 447 MPFR_ZIV_FREE (loop); 448 mpfr_clear (t); 449 mpfr_clear (q); 450 451 MPFR_SAVE_EXPO_FREE (expo); 452 return mpfr_check_range (y, res, rnd_mode); 453} 454 455int 456mpfr_cmp_q (mpfr_srcptr x, mpq_srcptr q) 457{ 458 mpfr_t t; 459 int res; 460 mpfr_prec_t p; 461 MPFR_SAVE_EXPO_DECL (expo); 462 463 /* GMP allows the user to set the denominator to 0. This is interpreted 464 by MPFR as the value being an infinity or NaN (probably better than 465 an assertion failure). */ 466 if (MPFR_UNLIKELY (mpz_sgn (mpq_denref (q)) == 0)) 467 { 468 /* q is an infinity or NaN */ 469 mpfr_flags_t old_flags; 470 471 mpfr_init2 (t, MPFR_PREC_MIN); 472 old_flags = __gmpfr_flags; 473 mpfr_set_q (t, q, MPFR_RNDN); 474 __gmpfr_flags = old_flags; 475 res = mpfr_cmp (x, t); 476 mpfr_clear (t); 477 return res; 478 } 479 480 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 481 return mpfr_cmp_si (x, mpq_sgn (q)); 482 483 MPFR_SAVE_EXPO_MARK (expo); 484 485 /* x < a/b ? <=> x*b < a */ 486 MPFR_MPZ_SIZEINBASE2 (p, mpq_denref (q)); 487 mpfr_init2 (t, MPFR_PREC(x) + p); 488 res = mpfr_mul_z (t, x, mpq_denref (q), MPFR_RNDN); 489 MPFR_ASSERTD (res == 0); 490 res = mpfr_cmp_z (t, mpq_numref (q)); 491 mpfr_clear (t); 492 493 MPFR_SAVE_EXPO_FREE (expo); 494 return res; 495} 496#endif 497 498#ifndef MPFR_USE_MINI_GMP 499int 500mpfr_cmp_f (mpfr_srcptr x, mpf_srcptr z) 501{ 502 mpfr_t t; 503 int res; 504 MPFR_SAVE_EXPO_DECL (expo); 505 506 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 507 return mpfr_cmp_si (x, mpf_sgn (z)); 508 509 MPFR_SAVE_EXPO_MARK (expo); 510 511 mpfr_init2 (t, MPFR_PREC_MIN + ABSIZ(z) * GMP_NUMB_BITS); 512 res = mpfr_set_f (t, z, MPFR_RNDN); 513 MPFR_ASSERTD (res == 0); 514 res = mpfr_cmp (x, t); 515 mpfr_clear (t); 516 517 MPFR_SAVE_EXPO_FREE (expo); 518 return res; 519} 520#endif 521