gmpy2_convert_mpfr.c
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 2 * gmpy_convert_mpfr.c * 3 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 4 * Python interface to the GMP or MPIR, MPFR, and MPC multiple precision * 5 * libraries. * 6 * * 7 * Copyright 2000 - 2009 Alex Martelli * 8 * * 9 * Copyright 2008 - 2021 Case Van Horsen * 10 * * 11 * This file is part of GMPY2. * 12 * * 13 * GMPY2 is free software: you can redistribute it and/or modify it under * 14 * the terms of the GNU Lesser General Public License as published by the * 15 * Free Software Foundation, either version 3 of the License, or (at your * 16 * option) any later version. * 17 * * 18 * GMPY2 is distributed in the hope that it will be useful, but WITHOUT * 19 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * 20 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public * 21 * License for more details. * 22 * * 23 * You should have received a copy of the GNU Lesser General Public * 24 * License along with GMPY2; if not, see <http://www.gnu.org/licenses/> * 25 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 26 27 /* This file contains all the conversion functions for MPFR data types. 28 * 29 * Overview 30 * -------- 31 * gmpy2 tries to optimize the performance and accuracy of conversions from 32 * other numeric types. The basic operations (+, -, *, /) are optimized to 33 * directly work with the basic types such as C longs or doubles. 34 * 35 * gmpy2 supports two different strategies for creating new references to 36 * an mpfr instance. If bits (or prec) is set to 0, the precison of the 37 * result exactly matches the precision of the context. The exponent range 38 * is also limited to the exponents defined in the contest. This is the 39 * default behavior of the mpfr() function. 40 * 41 * If bits (or prec) is set to 1, the precision of the result depends on the 42 * type of the source. 43 * 44 * If the source number is already a radix-2 floating point number, 45 * the precision is not changed. In practical terms, this only applies 46 * to sources operands that are either an mpfr or Python double. 47 * 48 * In addition, the exponent range is taken from the global emin/emax values 49 * set in the MPFR library. 50 * 51 */ 52 53 /* If prec == 0: 54 * Return an reference to a new mpfr instance. The context argument 55 * specifies the precision, exponent range, and whether or not subnormals 56 * are allowed. 57 * 58 * This is the default behavior of the mpfr() constructor. 59 * 60 * If prec == 1: 61 * Return an additional reference to an existing mpfr. The precision is not 62 * changed. The exponent is not checked and may be outside the bounds 63 * specified in the context argument. 64 * 65 * This is used internally for parsing arguments to functions. 66 * 67 * If prec >= 2: 68 * Return a reference to a new mpfr instance. The precision is taken from 69 * the argument. n instance with the specified precision and the 70 * exponent is valid in the current context. 71 * 72 * This is used by mpfr() with an optional precision keyword. 73 * 74 * Since references to existing objects may be returned, the result should not 75 * modified in-place. 76 */ 77 78 static MPFR_Object * 79 GMPy_MPFR_From_MPFR(MPFR_Object *obj, mpfr_prec_t prec, CTXT_Object *context) 80 { 81 MPFR_Object *result = NULL; 82 83 /* Optimize the critical case when prec==1 or obj is NaN or Inf. */ 84 85 if (prec == 1 || !mpfr_number_p(obj->f)) { 86 Py_INCREF((PyObject*)obj); 87 return obj; 88 } 89 90 CHECK_CONTEXT(context); 91 92 if (prec == 0) 93 prec = GET_MPFR_PREC(context); 94 95 /* Try to identify when an additional reference to existing instance can 96 * be returned. It is possible when (1) the precision matches, (2) the 97 * exponent is valid and not in the range that might require subnormal- 98 * ization, and (3) subnormalize is not enabled. 99 */ 100 101 if (prec == mpfr_get_prec(obj->f) && 102 !context->ctx.subnormalize && 103 obj->f->_mpfr_exp >= (context->ctx.emin + mpfr_get_prec(obj->f) - 1) && 104 obj->f->_mpfr_exp <= context->ctx.emax 105 ) { 106 107 Py_INCREF((PyObject*)obj); 108 return obj; 109 } 110 111 if ((result = GMPy_MPFR_New(prec, context))) { 112 mpfr_clear_flags(); 113 result->rc = mpfr_set(result->f, obj->f, GET_MPFR_ROUND(context)); 114 _GMPy_MPFR_Cleanup(&result, context); 115 } 116 return result; 117 } 118 119 static MPFR_Object * 120 GMPy_MPFR_From_PyIntOrLong(PyObject *obj, mpfr_prec_t prec, CTXT_Object *context) 121 { 122 MPFR_Object *result = NULL; 123 MPZ_Object *tempx = NULL; 124 int error, was_one = 0; 125 long temp; 126 127 CHECK_CONTEXT(context); 128 129 if (prec == 0) 130 prec = GET_MPFR_PREC(context); 131 132 if (prec == 1) { 133 /* This should be made accurate, but is it worth the overhead? */ 134 /* It is only used if the value fits in a C long. */ 135 prec = 64; 136 was_one = 1; 137 } 138 139 temp = GMPy_Integer_AsLongAndError(obj, &error); 140 141 if (error) { 142 if (!(tempx = GMPy_MPZ_From_PyIntOrLong(obj, context))) { 143 return NULL; 144 } 145 if (was_one) 146 prec = 1; 147 result = GMPy_MPFR_From_MPZ(tempx, prec, context); 148 Py_DECREF((PyObject*)tempx); 149 return result; 150 } 151 else { 152 if (!(result = GMPy_MPFR_New(prec, context))) 153 return NULL; 154 mpfr_clear_flags(); 155 result->rc = mpfr_set_si(result->f, temp, GET_MPFR_ROUND(context)); 156 if (!was_one) { 157 GMPY_MPFR_CHECK_RANGE(result, context); 158 } 159 GMPY_MPFR_EXCEPTIONS(result, context); 160 } 161 162 return result; 163 } 164 165 /* If prec==0, then the precision of the current context is used. 166 * 167 * If prec==1, then the precision of Python float type is used (typically 53). 168 * 169 * If prec>=2, then the specified precision is used. 170 */ 171 172 static MPFR_Object * 173 GMPy_MPFR_From_PyFloat(PyObject *obj, mpfr_prec_t prec, CTXT_Object *context) 174 { 175 MPFR_Object *result; 176 177 CHECK_CONTEXT(context); 178 179 if (prec == 0) 180 prec = GET_MPFR_PREC(context); 181 else if (prec == 1) 182 prec = DBL_MANT_DIG; 183 184 if ((result = GMPy_MPFR_New(prec, context))) { 185 mpfr_clear_flags(); 186 result->rc = mpfr_set_d(result->f, PyFloat_AS_DOUBLE(obj), GET_MPFR_ROUND(context)); 187 if (prec != 1) { 188 GMPY_MPFR_CHECK_RANGE(result, context); 189 } 190 GMPY_MPFR_SUBNORMALIZE(result, context); 191 GMPY_MPFR_EXCEPTIONS(result, context); 192 } 193 return result; 194 } 195 196 /* If prec==0, then the precision of the current context is used. 197 * 198 * If prec==1, then an exact conversion is done by using the bit length of the 199 * argument as the precision. 200 * 201 * If prec>=2, then the specified precision is used. 202 */ 203 204 static MPFR_Object * 205 GMPy_MPFR_From_MPZ(MPZ_Object *obj, mpfr_prec_t prec, CTXT_Object *context) 206 { 207 MPFR_Object *result; 208 int was_one = 0; 209 size_t bitlen; 210 211 CHECK_CONTEXT(context); 212 213 if (prec == 0) 214 prec = GET_MPFR_PREC(context); 215 216 if (prec == 1) { 217 bitlen = mpz_sizeinbase(obj->z, 2); 218 if (bitlen < MPFR_PREC_MIN) { 219 bitlen = MPFR_PREC_MIN; 220 } 221 if (bitlen > MPFR_PREC_MAX) { 222 OVERFLOW_ERROR("'mpz' to large to convert to 'mpfr'\n"); 223 return NULL; 224 } 225 prec = (mpfr_prec_t)bitlen; 226 was_one = 1; 227 } 228 229 if ((result = GMPy_MPFR_New(prec, context))) { 230 mpfr_clear_flags(); 231 result->rc = mpfr_set_z(result->f, obj->z, GET_MPFR_ROUND(context)); 232 if (!was_one) { 233 GMPY_MPFR_CHECK_RANGE(result, context); 234 } 235 GMPY_MPFR_EXCEPTIONS(result, context); 236 } 237 238 return result; 239 } 240 241 /* If prec<2, then the precision of the current context is used. 242 * 243 * If prec>=2, then the specified precision is used. 244 */ 245 246 static MPFR_Object * 247 GMPy_MPFR_From_MPQ(MPQ_Object *obj, mpfr_prec_t prec, CTXT_Object *context) 248 { 249 MPFR_Object *result; 250 251 CHECK_CONTEXT(context); 252 253 if (prec < 2) 254 prec = GET_MPFR_PREC(context); 255 256 if ((result = GMPy_MPFR_New(prec, context))) { 257 mpfr_clear_flags(); 258 result->rc = mpfr_set_q(result->f, obj->q, GET_MPFR_ROUND(context)); 259 GMPY_MPFR_CHECK_RANGE(result, context); 260 GMPY_MPFR_SUBNORMALIZE(result, context); 261 GMPY_MPFR_EXCEPTIONS(result, context); 262 } 263 264 return result; 265 } 266 267 /* If prec<2, then the precision of the current context is used. 268 * 269 * If prec>=2, then the specified precision is used. 270 */ 271 272 static MPFR_Object * 273 GMPy_MPFR_From_Fraction(PyObject *obj, mpfr_prec_t prec, CTXT_Object *context) 274 { 275 MPFR_Object *result = NULL; 276 MPQ_Object *tempq; 277 278 CHECK_CONTEXT(context); 279 280 if ((tempq = GMPy_MPQ_From_Fraction(obj, context))) { 281 result = GMPy_MPFR_From_MPQ(tempq, prec, context); 282 Py_DECREF((PyObject*)tempq); 283 } 284 285 return result; 286 } 287 288 static MPFR_Object * 289 GMPy_MPFR_From_PyStr(PyObject *s, int base, mpfr_prec_t prec, CTXT_Object *context) 290 { 291 MPFR_Object *result; 292 MPQ_Object *tempq; 293 char *cp, *endptr; 294 Py_ssize_t len; 295 PyObject *ascii_str = NULL; 296 297 CHECK_CONTEXT(context); 298 299 if (prec < 2) 300 prec = GET_MPFR_PREC(context); 301 302 if (PyBytes_Check(s)) { 303 len = PyBytes_Size(s); 304 cp = PyBytes_AsString(s); 305 } 306 else if (PyUnicode_Check(s)) { 307 ascii_str = PyUnicode_AsASCIIString(s); 308 if (!ascii_str) { 309 VALUE_ERROR("string contains non-ASCII characters"); 310 return NULL; 311 } 312 len = PyBytes_Size(ascii_str); 313 cp = PyBytes_AsString(ascii_str); 314 } 315 else { 316 TYPE_ERROR("object is not string or Unicode"); 317 return NULL; 318 } 319 320 /* Check for leading base indicators. */ 321 if (base == 0) { 322 if (len > 2 && cp[0] == '0') { 323 if (cp[1] == 'b') { base = 2; cp += 2; len -= 2; } 324 else if (cp[1] == 'x') { base = 16; cp += 2; len -= 2; } 325 else { base = 10; } 326 } 327 else { 328 base = 10; 329 } 330 } 331 else if (cp[0] == '0') { 332 /* If the specified base matches the leading base indicators, then 333 * we need to skip the base indicators. 334 */ 335 if (cp[1] =='b' && base == 2) { cp += 2; len -= 2; } 336 else if (cp[1] =='x' && base == 16) { cp += 2; len -= 2; } 337 } 338 339 340 if (!(result = GMPy_MPFR_New(prec, context))) { 341 Py_XDECREF(ascii_str); 342 return NULL; 343 } 344 345 /* delegate the rest to MPFR */ 346 mpfr_clear_flags(); 347 result->rc = mpfr_strtofr(result->f, cp, &endptr, base, GET_MPFR_ROUND(context)); 348 Py_XDECREF(ascii_str); 349 350 if (len != (Py_ssize_t)(endptr - cp)) { 351 VALUE_ERROR("invalid digits"); 352 Py_DECREF((PyObject*)result); 353 return NULL; 354 } 355 356 /* If the context requests subnormals and the result is in the range for subnormals, 357 * we use exact conversion via conversion to an mpq. 358 * 359 * The sticky bit returned by MPFR's string conversion appears to only reflect the 360 * portion of the string needed to compute the correctly rounded result. It does not 361 * accurately reflect whether or not the result is larger or smaller than the entire 362 * input string. A correct sticky bit is needed by mfpr_subnormalize. Converting the 363 * string to an mpq and then converting the mpq to an mpfr does properly set the 364 * sticky bit. 365 */ 366 367 if (base == 10 && 368 context->ctx.subnormalize && 369 result->f->_mpfr_exp <= context->ctx.emin + mpfr_get_prec(result->f) - 1) 370 { 371 372 if (!(tempq = GMPy_MPQ_From_PyStr(s, base, context))) { 373 Py_DECREF((PyObject*)result); 374 return NULL; 375 } 376 377 mpfr_clear_flags(); 378 result->rc = mpfr_set_q(result->f, tempq->q, GET_MPFR_ROUND(context)); 379 Py_DECREF((PyObject*)tempq); 380 } 381 382 GMPY_MPFR_CHECK_RANGE(result, context); 383 GMPY_MPFR_SUBNORMALIZE(result, context); 384 GMPY_MPFR_EXCEPTIONS(result, context); 385 386 return result; 387 } 388 389 /* GMPy_MPFR_From_Real() converts a real number (see IS_REAL()) to an mpfr. 390 * 391 * If prec==0, then the result has the precision of the current context. 392 * 393 * If prec==1 and the value can be converted exactly (i.e. the input value is 394 * a floating-point number using radix-2 representation or an integer), then 395 * the conversion is done with the maximum possible precision. If the input 396 * value can't be converted exactly, then the context precision is used. 397 * 398 * If prec >=2, then the specified precision is used. 399 * 400 * The return value is guaranteed to have a valid exponent. 401 */ 402 403 static MPFR_Object * 404 GMPy_MPFR_From_RealWithType(PyObject *obj, int xtype, mp_prec_t prec, CTXT_Object *context) 405 { 406 CHECK_CONTEXT(context); 407 408 if (IS_TYPE_MPFR(xtype)) 409 return GMPy_MPFR_From_MPFR((MPFR_Object*)obj, prec, context); 410 411 if (IS_TYPE_PyFloat(xtype)) 412 return GMPy_MPFR_From_PyFloat(obj, prec, context); 413 414 if (IS_TYPE_MPQ(xtype)) 415 return GMPy_MPFR_From_MPQ((MPQ_Object*)obj, prec, context); 416 417 if (IS_TYPE_MPZANY(xtype)) 418 return GMPy_MPFR_From_MPZ((MPZ_Object*)obj, prec, context); 419 420 if (IS_TYPE_PyInteger(xtype)) 421 return GMPy_MPFR_From_PyIntOrLong(obj, prec, context); 422 423 if (IS_TYPE_PyFraction(xtype)) 424 return GMPy_MPFR_From_Fraction(obj, prec, context); 425 426 if (IS_TYPE_HAS_MPFR(xtype)) { 427 MPFR_Object *res = (MPFR_Object *) PyObject_CallMethod(obj, "__mpfr__", NULL); 428 429 if (res != NULL && MPFR_Check(res)) { 430 return res; 431 } 432 else { 433 Py_XDECREF((PyObject*)res); 434 goto error; 435 } 436 } 437 438 if (IS_TYPE_HAS_MPQ(xtype)) { 439 MPQ_Object *res = (MPQ_Object *) PyObject_CallMethod(obj, "__mpq__", NULL); 440 441 if (res != NULL && MPQ_Check(res)) { 442 MPFR_Object * temp = GMPy_MPFR_From_MPQ(res, prec, context); 443 Py_DECREF(res); 444 return temp; 445 } 446 else { 447 Py_XDECREF((PyObject*)res); 448 goto error; 449 } 450 } 451 452 if (IS_TYPE_HAS_MPZ(xtype)) { 453 MPZ_Object *res = (MPZ_Object *) PyObject_CallMethod(obj, "__mpz__", NULL); 454 455 if (res != NULL && MPZ_Check(res)) { 456 MPFR_Object * temp = GMPy_MPFR_From_MPZ(res, prec, context); 457 Py_DECREF(res); 458 return temp; 459 } 460 else { 461 Py_XDECREF((PyObject*)res); 462 goto error; 463 } 464 } 465 466 error: 467 TYPE_ERROR("object could not be converted to 'mpfr'"); 468 return NULL; 469 } 470 471 static MPFR_Object * 472 GMPy_MPFR_From_Real(PyObject *obj, mp_prec_t prec, CTXT_Object *context) 473 { 474 return GMPy_MPFR_From_RealWithType(obj, GMPy_ObjectType(obj), 475 prec, context); 476 } 477 478 static MPFR_Object * 479 GMPy_MPFR_From_RealWithTypeAndCopy(PyObject *obj, int xtype, mp_prec_t prec, CTXT_Object *context) 480 { 481 MPFR_Object *result = NULL, *temp = NULL; 482 483 result = GMPy_MPFR_From_RealWithType(obj, xtype, prec, context); 484 485 if (result == NULL) 486 return result; 487 488 if (Py_REFCNT(result) == 1) 489 return result; 490 491 if (!(temp = GMPy_MPFR_New(mpfr_get_prec(result->f), context))) { 492 /* LCOV_EXCL_START */ 493 return NULL; 494 /* LCOV_EXCL_STOP */ 495 } 496 497 /* Since the precision of temp is the same as the precision of result, 498 * there shouldn't be any rounding. 499 */ 500 501 mpfr_set(temp->f, result->f, MPFR_RNDN); 502 Py_DECREF((PyObject*)result); 503 return temp; 504 } 505 506 static MPZ_Object * 507 GMPy_MPZ_From_MPFR(MPFR_Object *obj, CTXT_Object *context) 508 { 509 MPZ_Object *result; 510 511 CHECK_CONTEXT(context); 512 513 if ((result = GMPy_MPZ_New(context))) { 514 if (mpfr_nan_p(MPFR(obj))) { 515 Py_DECREF((PyObject*)result); 516 VALUE_ERROR("'mpz' does not support NaN"); 517 return NULL; 518 } 519 if (mpfr_inf_p(MPFR(obj))) { 520 Py_DECREF((PyObject*)result); 521 OVERFLOW_ERROR("'mpz' does not support Infinity"); 522 return NULL; 523 } 524 /* return code is ignored */ 525 mpfr_get_z(result->z, MPFR(obj), GET_MPFR_ROUND(context)); 526 } 527 528 return result; 529 } 530 531 static XMPZ_Object * 532 GMPy_XMPZ_From_MPFR(MPFR_Object *self, CTXT_Object *context) 533 { 534 XMPZ_Object *result; 535 536 CHECK_CONTEXT(context); 537 538 if ((result = GMPy_XMPZ_New(context))) { 539 if (mpfr_nan_p(MPFR(self))) { 540 Py_DECREF((PyObject*)result); 541 VALUE_ERROR("'xmpz' does not support NaN"); 542 return NULL; 543 } 544 if (mpfr_inf_p(MPFR(self))) { 545 Py_DECREF((PyObject*)result); 546 OVERFLOW_ERROR("'xmpz' does not support Infinity"); 547 return NULL; 548 } 549 /* return code is ignored */ 550 mpfr_get_z(result->z, MPFR(self), GET_MPFR_ROUND(context)); 551 } 552 553 return result; 554 } 555 556 /* Return the simpliest rational number that approximates 'self' to the 557 * requested error bound. 558 * 559 * 'bits' is used as the working precision for the calculations. 560 * - If (bits == 0) then the precision of 'self' is used. 561 * - If (bits > precision of 'self') then use precision of 'self'. 562 * 563 * 'err' is the maximum error in the rational approximation. 564 * - If (err > 0) then prec 565 * - If (err == NULL), then the requested precision is 1/(2**prec). 566 * This should return the smallest fraction that returns the 567 * same interval that includes 'self'. 568 * 'err'. If 'err' is negative, then the requested 569 * precision is -2**abs(int(err)). If 'err' is NULL, then the requested 570 * precision is -2**prec. If 'prec' is 0, then the requested precision is 571 * the precision of 'self'. 572 */ 573 574 575 static PyObject * 576 stern_brocot(MPFR_Object* self, MPFR_Object *err, mpfr_prec_t bits, int mayz, CTXT_Object *context) 577 { 578 PyObject* result = NULL; 579 MPQ_Object *mpqres = NULL; 580 MPZ_Object *mpzres = NULL; 581 582 int i, negative, errsign; 583 mpfr_t f, al, a, r1[3], r2[3], temperr, minerr, curerr, newerr, temp; 584 585 if (mpfr_nan_p(self->f)) { 586 VALUE_ERROR("Cannot convert NaN to a number."); 587 return NULL; 588 } 589 590 if (mpfr_inf_p(self->f)) { 591 OVERFLOW_ERROR("Cannot convert Infinity to a number."); 592 return NULL; 593 } 594 595 if (!bits) { 596 bits = mpfr_get_prec(self->f); 597 } 598 599 if (err) { 600 /* Make a copy. */ 601 mpfr_init2(temperr, mpfr_get_prec(err->f)); 602 mpfr_set(temperr, err->f, MPFR_RNDN); 603 } 604 else { 605 mpfr_init2(temperr, mpfr_get_prec(self->f)); 606 mpfr_set_ui(temperr, 0, MPFR_RNDN); 607 } 608 609 errsign = mpfr_sgn(temperr); 610 if (errsign <= 0 && (bits < 2 || bits > mpfr_get_prec(self->f))) { 611 VALUE_ERROR("Requested precision out-of-bounds."); 612 mpfr_clear(temperr); 613 return NULL; 614 } 615 616 if (errsign == 0) { 617 mpfr_set_si(temperr, 1, MPFR_RNDN); 618 mpfr_div_2exp(temperr, temperr, bits, MPFR_RNDN); 619 } 620 else if (errsign < 0) { 621 long ubits; 622 mpfr_abs(temperr, temperr, MPFR_RNDN); 623 mpfr_floor(temperr, temperr); 624 ubits = mpfr_get_si(temperr, MPFR_RNDN); 625 if (ubits < 2 || ubits > mpfr_get_prec(self->f)) { 626 VALUE_ERROR("Requested precision out-of-bounds."); 627 mpfr_clear(temperr); 628 return NULL; 629 } 630 mpfr_set_si(temperr, 1, MPFR_RNDN); 631 mpfr_div_2exp(temperr, temperr, ubits, MPFR_RNDN); 632 } 633 634 if (!(mpqres = GMPy_MPQ_New(context)) || 635 !(mpzres = GMPy_MPZ_New(context))) { 636 Py_XDECREF((PyObject*)mpqres); 637 Py_XDECREF((PyObject*)mpzres); 638 mpfr_clear(temperr); 639 return NULL; 640 } 641 642 mpfr_init2(minerr, bits); 643 mpfr_set(minerr, temperr, MPFR_RNDN); 644 mpfr_clear(temperr); 645 646 mpfr_init2(f, bits); 647 if (mpfr_sgn(self->f) < 0) { 648 negative = 1; 649 mpfr_abs(f, self->f, MPFR_RNDN); 650 } 651 else { 652 negative = 0; 653 mpfr_set(f, self->f, MPFR_RNDN); 654 } 655 656 mpfr_init2(al, bits); 657 mpfr_set(al, f, MPFR_RNDN); 658 mpfr_init2(a, bits); 659 mpfr_floor(a, al); 660 mpfr_init2(temp, bits); 661 for (i=0; i<3; ++i) { 662 mpfr_init2(r1[i], bits); 663 mpfr_init2(r2[i], bits); 664 } 665 666 mpfr_set_si(r1[0], 0, MPFR_RNDN); 667 mpfr_set_si(r1[1], 0, MPFR_RNDN); 668 mpfr_set_si(r1[2], 1, MPFR_RNDN); 669 mpfr_set_si(r2[0], 0, MPFR_RNDN); 670 mpfr_set_si(r2[1], 1, MPFR_RNDN); 671 mpfr_set(r2[2], a, MPFR_RNDN); 672 mpfr_init2(curerr, bits); 673 mpfr_init2(newerr, bits); 674 mpfr_reldiff(curerr, f, a, MPFR_RNDN); 675 676 while (mpfr_cmp(curerr, minerr) > 0) { 677 mpfr_sub(temp, al, a, MPFR_RNDN); 678 mpfr_ui_div(al, 1, temp, MPFR_RNDN); 679 mpfr_floor(a, al); 680 mpfr_swap(r1[0], r1[1]); 681 mpfr_swap(r1[1], r1[2]); 682 //mpfr_fma(r1[2], r1[1], a, r1[0], MPFR_RNDN); 683 mpfr_mul(r1[2], r1[1], a, MPFR_RNDN); 684 mpfr_add(r1[2], r1[2], r1[0], MPFR_RNDN); 685 mpfr_swap(r2[0], r2[1]); 686 mpfr_swap(r2[1], r2[2]); 687 //mpfr_fma(r2[2], r2[1], a, r2[0], MPFR_RNDN); 688 mpfr_mul(r2[2], r2[1], a, MPFR_RNDN); 689 mpfr_add(r2[2], r2[2], r2[0], MPFR_RNDN); 690 mpfr_div(temp, r2[2], r1[2], MPFR_RNDN); 691 mpfr_reldiff(newerr, f, temp, MPFR_RNDN); 692 if(mpfr_cmp(curerr, newerr) <= 0) { 693 mpfr_swap(r1[1],r1[2]); 694 mpfr_swap(r2[1],r2[2]); 695 break; 696 } 697 mpfr_swap(curerr, newerr); 698 } 699 700 /* Note: both mpqres and mpzrec have been created. Remember to delete the 701 * one you don't need. 702 */ 703 704 if (mayz && (mpfr_cmp_ui(r1[2],1) == 0)) { 705 mpfr_get_z(mpzres->z, r2[2], MPFR_RNDN); 706 if (negative) { 707 mpz_neg(mpzres->z, mpzres->z); 708 } 709 result = (PyObject*)mpzres; 710 Py_DECREF(mpqres); 711 } 712 else { 713 mpfr_get_z(mpq_numref(mpqres->q), r2[2], MPFR_RNDN); 714 mpfr_get_z(mpq_denref(mpqres->q), r1[2], MPFR_RNDN); 715 if (negative) { 716 mpz_neg(mpq_numref(mpqres->q), mpq_numref(mpqres->q)); 717 } 718 result = (PyObject*)mpqres; 719 Py_DECREF(mpzres); 720 } 721 722 mpfr_clear(minerr); 723 mpfr_clear(al); 724 mpfr_clear(a); 725 mpfr_clear(f); 726 for (i=0; i<3; ++i) { 727 mpfr_clear(r1[i]); 728 mpfr_clear(r2[i]); 729 } 730 mpfr_clear(curerr); 731 mpfr_clear(newerr); 732 mpfr_clear(temp); 733 return result; 734 } 735 736 static MPQ_Object * 737 GMPy_MPQ_From_MPFR(MPFR_Object *self, CTXT_Object *context) 738 { 739 mpfr_exp_t temp, twocount; 740 MPQ_Object *result; 741 742 CHECK_CONTEXT(context); 743 744 if (mpfr_nan_p(self->f)) { 745 VALUE_ERROR("can not convert NaN to MPQ"); 746 return NULL; 747 } 748 749 if (mpfr_inf_p(self->f)) { 750 OVERFLOW_ERROR("can not convert Infinity to MPQ"); 751 return NULL; 752 } 753 754 if (!(result = GMPy_MPQ_New(context))) { 755 return NULL; 756 } 757 758 if (mpfr_zero_p(self->f)) { 759 mpz_set_ui(mpq_numref(result->q), 0); 760 mpz_set_ui(mpq_denref(result->q), 1); 761 } 762 else { 763 temp = mpfr_get_z_2exp(mpq_numref(result->q), self->f); 764 twocount = (mpfr_exp_t)mpz_scan1(mpq_numref(result->q), 0); 765 if (twocount) { 766 temp += twocount; 767 mpz_div_2exp(mpq_numref(result->q), mpq_numref(result->q), twocount); 768 } 769 mpz_set_ui(mpq_denref(result->q), 1); 770 if (temp > 0) 771 mpz_mul_2exp(mpq_numref(result->q), mpq_numref(result->q), temp); 772 else if (temp < 0) 773 mpz_mul_2exp(mpq_denref(result->q), mpq_denref(result->q), -temp); 774 } 775 return result; 776 } 777 778 static PyObject * 779 GMPy_PyIntOrLong_From_MPFR(MPFR_Object *obj, CTXT_Object *context) 780 { 781 PyObject *result; 782 MPZ_Object *tempz; 783 784 CHECK_CONTEXT(context); 785 786 if (!(tempz = GMPy_MPZ_From_MPFR(obj, context))) 787 return NULL; 788 789 result = GMPy_PyIntOrLong_From_MPZ(tempz, context); 790 Py_DECREF((PyObject*)tempz); 791 792 return result; 793 } 794 795 static PyObject * 796 GMPy_MPFR_Int_Slot(MPFR_Object *self) { 797 return GMPy_PyIntOrLong_From_MPFR(self, NULL); 798 } 799 800 #ifdef PY2 801 static PyObject * 802 GMPy_PyLong_From_MPFR(MPFR_Object *obj, CTXT_Object *context) 803 { 804 PyObject *result; 805 MPZ_Object *tempz; 806 807 CHECK_CONTEXT(context); 808 809 if (!(tempz = GMPy_MPZ_From_MPFR(obj, context))) 810 return NULL; 811 812 result = GMPy_PyLong_From_MPZ(tempz, context); 813 Py_DECREF((PyObject*)tempz); 814 815 return result; 816 } 817 818 static PyObject * 819 GMPy_MPFR_Long_Slot(MPFR_Object *self) { 820 return GMPy_PyLong_From_MPFR(self, NULL); 821 } 822 #endif 823 824 static PyObject * 825 GMPy_PyFloat_From_MPFR(MPFR_Object *self, CTXT_Object *context) 826 { 827 double res; 828 829 CHECK_CONTEXT(context); 830 831 res = mpfr_get_d(self->f, GET_MPFR_ROUND(context)); 832 833 return PyFloat_FromDouble(res); 834 } 835 836 static PyObject * 837 GMPy_MPFR_Float_Slot(MPFR_Object *self) { 838 return GMPy_PyFloat_From_MPFR(self, NULL); 839 } 840 841 static PyObject* 842 GMPy_PyStr_From_MPFR(MPFR_Object *self, int base, int digits, CTXT_Object *context) 843 { 844 PyObject *result; 845 char *buffer; 846 mpfr_exp_t the_exp; 847 848 CHECK_CONTEXT(context); 849 850 /* check arguments are valid */ 851 if (!((base >= 2) && (base <= 62))) { 852 VALUE_ERROR("base must be in the interval [2,62]"); 853 return NULL; 854 } 855 if ((digits < 0) || (digits == 1)) { 856 VALUE_ERROR("digits must be 0 or >= 2"); 857 return NULL; 858 } 859 860 /* Process special cases first */ 861 if (!(mpfr_regular_p(self->f))) { 862 if (mpfr_nan_p(self->f)) { 863 return Py_BuildValue("(sii)", "nan", 0, 0); 864 } 865 else if (mpfr_inf_p(self->f) && !mpfr_signbit(self->f)) { 866 return Py_BuildValue("(sii)", "inf", 0, 0); 867 } 868 else if (mpfr_inf_p(self->f) && mpfr_signbit(self->f)) { 869 return Py_BuildValue("(sii)", "-inf", 0, 0); 870 } 871 /* 0 is not considered a 'regular" number */ 872 else if (mpfr_signbit(self->f)) { 873 return Py_BuildValue("(sii)", "-0", 0, mpfr_get_prec(self->f)); 874 } 875 else { 876 return Py_BuildValue("(sii)", "0", 0, mpfr_get_prec(self->f)); 877 } 878 } 879 880 /* obtain digits-string and exponent */ 881 buffer = mpfr_get_str(0, &the_exp, base, digits, self->f, GET_MPFR_ROUND(context)); 882 if (!*buffer) { 883 SYSTEM_ERROR("Internal error in Pympfr_To_PyStr"); 884 return NULL; 885 } 886 887 result = Py_BuildValue("(sii)", buffer, the_exp, mpfr_get_prec(self->f)); 888 mpfr_free_str(buffer); 889 return result; 890 } 891 892 #ifdef SHARED 893 /* 894 * coerce any number to a mpf 895 */ 896 897 int 898 GMPy_MPFR_ConvertArg(PyObject *arg, PyObject **ptr) 899 { 900 MPFR_Object* newob = GMPy_MPFR_From_RealWithType(arg, GMPy_ObjectType(arg), 901 1, NULL); 902 903 if (newob) { 904 *ptr = (PyObject*)newob; 905 return 1; 906 } 907 else { 908 TYPE_ERROR("argument can not be converted to 'mpfr'"); 909 return 0; 910 } 911 } 912 #endif 913 914 /* str and repr implementations for mpfr */ 915 static PyObject * 916 GMPy_MPFR_Str_Slot(MPFR_Object *self) 917 { 918 PyObject *result, *temp; 919 long precision; 920 char fmtstr[60]; 921 922 precision = (long)(log10(2) * (double)mpfr_get_prec(MPFR(self))) + 2; 923 924 sprintf(fmtstr, "{0:.%ldg}", precision); 925 926 temp = Py_BuildValue("s", fmtstr); 927 if (!temp) 928 return NULL; 929 result = PyObject_CallMethod(temp, "format", "O", self); 930 Py_DECREF(temp); 931 return result; 932 } 933 934 static PyObject * 935 GMPy_MPFR_Repr_Slot(MPFR_Object *self) 936 { 937 PyObject *result, *temp; 938 long precision, bits; 939 char fmtstr[60]; 940 941 bits = mpfr_get_prec(MPFR(self)); 942 precision = (long)(log10(2) * (double)bits) + 2; 943 944 if (mpfr_number_p(MPFR(self)) && bits != DBL_MANT_DIG) 945 sprintf(fmtstr, "mpfr('{0:.%ldg}',%ld)", precision, bits); 946 else 947 sprintf(fmtstr, "mpfr('{0:.%ldg}')", precision); 948 949 temp = Py_BuildValue("s", fmtstr); 950 if (!temp) 951 return NULL; 952 result = PyObject_CallMethod(temp, "format", "O", self); 953 Py_DECREF(temp); 954 return result; 955 } 956 957 static PyObject * 958 mpfr_ascii(mpfr_t self, int base, int digits, int round) 959 { 960 PyObject *result; 961 char *buffer; 962 mpfr_exp_t the_exp; 963 964 /* Process special cases first */ 965 if (!(mpfr_regular_p(self))) { 966 if (mpfr_nan_p(self)) { 967 return Py_BuildValue("(sii)", "nan", 0, 0); 968 } 969 else if (mpfr_inf_p(self) && !mpfr_signbit(self)) { 970 return Py_BuildValue("(sii)", "inf", 0, 0); 971 } 972 else if (mpfr_inf_p(self) && mpfr_signbit(self)) { 973 return Py_BuildValue("(sii)", "-inf", 0, 0); 974 } 975 /* 0 is not considered a 'regular" number */ 976 else if (mpfr_signbit(self)) { 977 return Py_BuildValue("(sii)", "-0", 0, mpfr_get_prec(self)); 978 } 979 else { 980 return Py_BuildValue("(sii)", "0", 0, mpfr_get_prec(self)); 981 } 982 } 983 984 /* obtain digits-string and exponent */ 985 buffer = mpfr_get_str(0, &the_exp, base, digits, self, round); 986 if (!*buffer) { 987 SYSTEM_ERROR("Internal error in mpfr_ascii"); 988 return NULL; 989 } 990 991 result = Py_BuildValue("(sii)", buffer, the_exp, mpfr_get_prec(self)); 992 mpfr_free_str(buffer); 993 return result; 994 }