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