imrat.c (27782B)
1 /* 2 ** imrat.c 3 ** Arbitrary precision rational arithmetic routines. 4 ** See Copyright Notice in klisp.h 5 */ 6 7 /* 8 ** SOURCE NOTE: This is mostly from the IMath library, written by 9 ** M.J. Fromberger. It is adapted to klisp, mainly in the use of the 10 ** klisp allocator and fixing of digit size to 32 bits. 11 ** Imported from version (1.15) updated 01-Feb-2011 at 03:10 PM. 12 */ 13 14 #include "imrat.h" 15 #include <stdlib.h> 16 #include <string.h> 17 #include <ctype.h> 18 #include <assert.h> 19 20 /* Andres Navarro: klisp includes */ 21 #include "kobject.h" 22 #include "kstate.h" 23 #include "kmem.h" 24 #include "kerror.h" 25 26 /* {{{ Useful macros */ 27 28 #define TEMP(K) (temp + (K)) 29 #define SETUP(E, C) \ 30 do{if((res = (E)) != MP_OK) goto CLEANUP; ++(C);}while(0) 31 32 /* Argument checking: 33 Use CHECK() where a return value is required; NRCHECK() elsewhere */ 34 #define CHECK(TEST) assert(TEST) 35 #define NRCHECK(TEST) assert(TEST) 36 37 /* }}} */ 38 39 /* Reduce the given rational, in place, to lowest terms and canonical 40 form. Zero is represented as 0/1, one as 1/1. Signs are adjusted 41 so that the sign of the numerator is definitive. */ 42 static mp_result s_rat_reduce(klisp_State *K, mp_rat r); 43 44 /* Common code for addition and subtraction operations on rationals. */ 45 static mp_result s_rat_combine(klisp_State *K, mp_rat a, mp_rat b, mp_rat c, 46 mp_result (*comb_f) 47 (klisp_State *,mp_int, mp_int, mp_int)); 48 49 /* {{{ mp_rat_init(r) */ 50 51 mp_result mp_rat_init(klisp_State *K, mp_rat r) 52 { 53 return mp_rat_init_size(K, r, 0, 0); 54 } 55 56 /* }}} */ 57 58 /* {{{ mp_rat_alloc() */ 59 60 mp_rat mp_rat_alloc(klisp_State *K) 61 { 62 mp_rat out = klispM_new(K, mpq_t); 63 64 if(out != NULL) { 65 if(mp_rat_init(K, out) != MP_OK) { 66 klispM_free(K, out); 67 return NULL; 68 } 69 } 70 71 return out; 72 } 73 74 /* }}} */ 75 76 /* {{{ mp_rat_init_size(r, n_prec, d_prec) */ 77 78 mp_result mp_rat_init_size(klisp_State *K, mp_rat r, mp_size n_prec, 79 mp_size d_prec) 80 { 81 mp_result res; 82 83 if((res = mp_int_init_size(K, MP_NUMER_P(r), n_prec)) != MP_OK) 84 return res; 85 if((res = mp_int_init_size(K, MP_DENOM_P(r), d_prec)) != MP_OK) { 86 mp_int_clear(K, MP_NUMER_P(r)); 87 return res; 88 } 89 90 return mp_int_set_value(K, MP_DENOM_P(r), 1); 91 } 92 93 /* }}} */ 94 95 /* {{{ mp_rat_init_copy(r, old) */ 96 97 mp_result mp_rat_init_copy(klisp_State *K, mp_rat r, mp_rat old) 98 { 99 mp_result res; 100 101 if((res = mp_int_init_copy(K, MP_NUMER_P(r), MP_NUMER_P(old))) != MP_OK) 102 return res; 103 if((res = mp_int_init_copy(K, MP_DENOM_P(r), MP_DENOM_P(old))) != MP_OK) 104 mp_int_clear(K, MP_NUMER_P(r)); 105 106 return res; 107 } 108 109 /* }}} */ 110 111 /* {{{ mp_rat_set_value(r, numer, denom) */ 112 113 mp_result mp_rat_set_value(klisp_State *K, mp_rat r, int numer, int denom) 114 { 115 mp_result res; 116 117 if(denom == 0) 118 return MP_UNDEF; 119 120 if((res = mp_int_set_value(K, MP_NUMER_P(r), numer)) != MP_OK) 121 return res; 122 if((res = mp_int_set_value(K, MP_DENOM_P(r), denom)) != MP_OK) 123 return res; 124 125 return s_rat_reduce(K, r); 126 } 127 128 /* }}} */ 129 130 /* {{{ mp_rat_clear(r) */ 131 132 void mp_rat_clear(klisp_State *K, mp_rat r) 133 { 134 mp_int_clear(K, MP_NUMER_P(r)); 135 mp_int_clear(K, MP_DENOM_P(r)); 136 } 137 138 /* }}} */ 139 140 /* {{{ mp_rat_free(r) */ 141 142 void mp_rat_free(klisp_State *K, mp_rat r) 143 { 144 NRCHECK(r != NULL); 145 146 if(r->num.digits != NULL) 147 mp_rat_clear(K, r); 148 149 klispM_free(K, r); 150 } 151 152 /* }}} */ 153 154 /* {{{ mp_rat_numer(r, z) */ 155 156 mp_result mp_rat_numer(klisp_State *K, mp_rat r, mp_int z) 157 { 158 return mp_int_copy(K, MP_NUMER_P(r), z); 159 } 160 161 /* }}} */ 162 163 /* {{{ mp_rat_denom(r, z) */ 164 165 mp_result mp_rat_denom(klisp_State *K, mp_rat r, mp_int z) 166 { 167 return mp_int_copy(K, MP_DENOM_P(r), z); 168 } 169 170 /* }}} */ 171 172 /* {{{ mp_rat_sign(r) */ 173 174 mp_sign mp_rat_sign(mp_rat r) 175 { 176 return MP_SIGN(MP_NUMER_P(r)); 177 } 178 179 /* }}} */ 180 181 /* {{{ mp_rat_copy(a, c) */ 182 183 mp_result mp_rat_copy(klisp_State *K, mp_rat a, mp_rat c) 184 { 185 mp_result res; 186 187 if((res = mp_int_copy(K, MP_NUMER_P(a), MP_NUMER_P(c))) != MP_OK) 188 return res; 189 190 res = mp_int_copy(K, MP_DENOM_P(a), MP_DENOM_P(c)); 191 return res; 192 } 193 194 /* }}} */ 195 196 /* {{{ mp_rat_zero(r) */ 197 198 void mp_rat_zero(klisp_State *K, mp_rat r) 199 { 200 mp_int_zero(MP_NUMER_P(r)); 201 mp_int_set_value(K, MP_DENOM_P(r), 1); 202 } 203 204 /* }}} */ 205 206 /* {{{ mp_rat_abs(a, c) */ 207 208 mp_result mp_rat_abs(klisp_State *K, mp_rat a, mp_rat c) 209 { 210 mp_result res; 211 212 if((res = mp_int_abs(K, MP_NUMER_P(a), MP_NUMER_P(c))) != MP_OK) 213 return res; 214 215 res = mp_int_abs(K, MP_DENOM_P(a), MP_DENOM_P(c)); 216 return res; 217 } 218 219 /* }}} */ 220 221 /* {{{ mp_rat_neg(a, c) */ 222 223 mp_result mp_rat_neg(klisp_State *K, mp_rat a, mp_rat c) 224 { 225 mp_result res; 226 227 if((res = mp_int_neg(K, MP_NUMER_P(a), 228 MP_NUMER_P(c))) != MP_OK) 229 return res; 230 231 res = mp_int_copy(K, MP_DENOM_P(a), MP_DENOM_P(c)); 232 return res; 233 } 234 235 /* }}} */ 236 237 /* {{{ mp_rat_recip(a, c) */ 238 239 mp_result mp_rat_recip(klisp_State *K, mp_rat a, mp_rat c) 240 { 241 mp_result res; 242 243 if(mp_rat_compare_zero(a) == 0) 244 return MP_UNDEF; 245 246 if((res = mp_rat_copy(K, a, c)) != MP_OK) 247 return res; 248 249 mp_int_swap(MP_NUMER_P(c), MP_DENOM_P(c)); 250 251 /* Restore the signs of the swapped elements */ 252 { 253 mp_sign tmp = MP_SIGN(MP_NUMER_P(c)); 254 255 MP_SIGN(MP_NUMER_P(c)) = MP_SIGN(MP_DENOM_P(c)); 256 MP_SIGN(MP_DENOM_P(c)) = tmp; 257 } 258 259 return MP_OK; 260 } 261 262 /* }}} */ 263 264 /* {{{ mp_rat_add(a, b, c) */ 265 266 mp_result mp_rat_add(klisp_State *K, mp_rat a, mp_rat b, mp_rat c) 267 { 268 return s_rat_combine(K, a, b, c, mp_int_add); 269 270 } 271 272 /* }}} */ 273 274 /* {{{ mp_rat_sub(a, b, c) */ 275 276 mp_result mp_rat_sub(klisp_State *K, mp_rat a, mp_rat b, mp_rat c) 277 { 278 return s_rat_combine(K, a, b, c, mp_int_sub); 279 280 } 281 282 /* }}} */ 283 284 /* {{{ mp_rat_mul(a, b, c) */ 285 286 mp_result mp_rat_mul(klisp_State *K, mp_rat a, mp_rat b, mp_rat c) 287 { 288 mp_result res; 289 290 if((res = mp_int_mul(K, MP_NUMER_P(a), MP_NUMER_P(b), 291 MP_NUMER_P(c))) != MP_OK) 292 return res; 293 294 if(mp_int_compare_zero(MP_NUMER_P(c)) != 0) { 295 if((res = mp_int_mul(K, MP_DENOM_P(a), MP_DENOM_P(b), 296 MP_DENOM_P(c))) != MP_OK) 297 return res; 298 } 299 300 return s_rat_reduce(K, c); 301 } 302 303 /* }}} */ 304 305 /* {{{ mp_int_div(a, b, c) */ 306 307 mp_result mp_rat_div(klisp_State *K, mp_rat a, mp_rat b, mp_rat c) 308 { 309 mp_result res = MP_OK; 310 311 if(mp_rat_compare_zero(b) == 0) 312 return MP_UNDEF; 313 314 if(c == a || c == b) { 315 mpz_t tmp; 316 317 if((res = mp_int_init(&tmp)) != MP_OK) return res; 318 if((res = mp_int_mul(K, MP_NUMER_P(a), MP_DENOM_P(b), &tmp)) != MP_OK) 319 goto CLEANUP; 320 if((res = mp_int_mul(K, MP_DENOM_P(a), MP_NUMER_P(b), 321 MP_DENOM_P(c))) != MP_OK) 322 goto CLEANUP; 323 res = mp_int_copy(K, &tmp, MP_NUMER_P(c)); 324 325 CLEANUP: 326 mp_int_clear(K, &tmp); 327 } 328 else { 329 if((res = mp_int_mul(K, MP_NUMER_P(a), MP_DENOM_P(b), 330 MP_NUMER_P(c))) != MP_OK) 331 return res; 332 if((res = mp_int_mul(K, MP_DENOM_P(a), MP_NUMER_P(b), 333 MP_DENOM_P(c))) != MP_OK) 334 return res; 335 } 336 337 if(res != MP_OK) 338 return res; 339 else 340 return s_rat_reduce(K, c); 341 } 342 343 /* }}} */ 344 345 /* {{{ mp_rat_add_int(a, b, c) */ 346 347 mp_result mp_rat_add_int(klisp_State *K, mp_rat a, mp_int b, mp_rat c) 348 { 349 mpz_t tmp; 350 mp_result res; 351 352 if((res = mp_int_init_copy(K, &tmp, b)) != MP_OK) 353 return res; 354 355 if((res = mp_int_mul(K, &tmp, MP_DENOM_P(a), &tmp)) != MP_OK) 356 goto CLEANUP; 357 358 if((res = mp_rat_copy(K, a, c)) != MP_OK) 359 goto CLEANUP; 360 361 if((res = mp_int_add(K, MP_NUMER_P(c), &tmp, MP_NUMER_P(c))) != MP_OK) 362 goto CLEANUP; 363 364 res = s_rat_reduce(K, c); 365 366 CLEANUP: 367 mp_int_clear(K, &tmp); 368 return res; 369 } 370 371 /* }}} */ 372 373 /* {{{ mp_rat_sub_int(a, b, c) */ 374 375 mp_result mp_rat_sub_int(klisp_State *K, mp_rat a, mp_int b, mp_rat c) 376 { 377 mpz_t tmp; 378 mp_result res; 379 380 if((res = mp_int_init_copy(K, &tmp, b)) != MP_OK) 381 return res; 382 383 if((res = mp_int_mul(K, &tmp, MP_DENOM_P(a), &tmp)) != MP_OK) 384 goto CLEANUP; 385 386 if((res = mp_rat_copy(K, a, c)) != MP_OK) 387 goto CLEANUP; 388 389 if((res = mp_int_sub(K, MP_NUMER_P(c), &tmp, MP_NUMER_P(c))) != MP_OK) 390 goto CLEANUP; 391 392 res = s_rat_reduce(K, c); 393 394 CLEANUP: 395 mp_int_clear(K, &tmp); 396 return res; 397 } 398 399 /* }}} */ 400 401 /* {{{ mp_rat_mul_int(a, b, c) */ 402 403 mp_result mp_rat_mul_int(klisp_State *K, mp_rat a, mp_int b, mp_rat c) 404 { 405 mp_result res; 406 407 if((res = mp_rat_copy(K, a, c)) != MP_OK) 408 return res; 409 410 if((res = mp_int_mul(K, MP_NUMER_P(c), b, MP_NUMER_P(c))) != MP_OK) 411 return res; 412 413 return s_rat_reduce(K, c); 414 } 415 416 /* }}} */ 417 418 /* {{{ mp_rat_div_int(a, b, c) */ 419 420 mp_result mp_rat_div_int(klisp_State *K, mp_rat a, mp_int b, mp_rat c) 421 { 422 mp_result res; 423 424 if(mp_int_compare_zero(b) == 0) 425 return MP_UNDEF; 426 427 if((res = mp_rat_copy(K, a, c)) != MP_OK) 428 return res; 429 430 if((res = mp_int_mul(K, MP_DENOM_P(c), b, MP_DENOM_P(c))) != MP_OK) 431 return res; 432 433 return s_rat_reduce(K, c); 434 } 435 436 /* }}} */ 437 438 /* {{{ mp_rat_expt(a, b, c) */ 439 440 mp_result mp_rat_expt(klisp_State *K, mp_rat a, mp_small b, mp_rat c) 441 { 442 mp_result res; 443 444 /* Special cases for easy powers. */ 445 if(b == 0) 446 return mp_rat_set_value(K, c, 1, 1); 447 else if(b == 1) 448 return mp_rat_copy(K, a, c); 449 450 /* Since rationals are always stored in lowest terms, it is not 451 necessary to reduce again when raising to an integer power. */ 452 if((res = mp_int_expt(K, MP_NUMER_P(a), b, MP_NUMER_P(c))) != MP_OK) 453 return res; 454 455 return mp_int_expt(K, MP_DENOM_P(a), b, MP_DENOM_P(c)); 456 } 457 458 /* }}} */ 459 460 /* {{{ mp_rat_compare(a, b) */ 461 462 int mp_rat_compare(klisp_State *K, mp_rat a, mp_rat b) 463 { 464 /* Quick check for opposite signs. Works because the sign of the 465 numerator is always definitive. */ 466 if(MP_SIGN(MP_NUMER_P(a)) != MP_SIGN(MP_NUMER_P(b))) { 467 if(MP_SIGN(MP_NUMER_P(a)) == MP_ZPOS) 468 return 1; 469 else 470 return -1; 471 } 472 else { 473 /* Compare absolute magnitudes; if both are positive, the answer 474 stands, otherwise it needs to be reflected about zero. */ 475 int cmp = mp_rat_compare_unsigned(K, a, b); 476 477 if(MP_SIGN(MP_NUMER_P(a)) == MP_ZPOS) 478 return cmp; 479 else 480 return -cmp; 481 } 482 } 483 484 /* }}} */ 485 486 /* {{{ mp_rat_compare_unsigned(a, b) */ 487 488 int mp_rat_compare_unsigned(klisp_State *K, mp_rat a, mp_rat b) 489 { 490 /* If the denominators are equal, we can quickly compare numerators 491 without multiplying. Otherwise, we actually have to do some work. */ 492 if(mp_int_compare_unsigned(MP_DENOM_P(a), MP_DENOM_P(b)) == 0) 493 return mp_int_compare_unsigned(MP_NUMER_P(a), MP_NUMER_P(b)); 494 495 else { 496 mpz_t temp[2]; 497 mp_result res; 498 int cmp = INT_MAX, last = 0; 499 500 /* t0 = num(a) * den(b), t1 = num(b) * den(a) */ 501 SETUP(mp_int_init_copy(K, TEMP(last), MP_NUMER_P(a)), last); 502 SETUP(mp_int_init_copy(K, TEMP(last), MP_NUMER_P(b)), last); 503 504 if((res = mp_int_mul(K, TEMP(0), MP_DENOM_P(b), TEMP(0))) != MP_OK || 505 (res = mp_int_mul(K, TEMP(1), MP_DENOM_P(a), TEMP(1))) != MP_OK) 506 goto CLEANUP; 507 508 cmp = mp_int_compare_unsigned(TEMP(0), TEMP(1)); 509 510 CLEANUP: 511 while(--last >= 0) 512 mp_int_clear(K, TEMP(last)); 513 514 return cmp; 515 } 516 } 517 518 /* }}} */ 519 520 /* {{{ mp_rat_compare_zero(r) */ 521 522 int mp_rat_compare_zero(mp_rat r) 523 { 524 return mp_int_compare_zero(MP_NUMER_P(r)); 525 526 } 527 528 /* }}} */ 529 530 /* {{{ mp_rat_compare_value(r, n, d) */ 531 532 int mp_rat_compare_value(klisp_State *K, mp_rat r, mp_small n, mp_small d) 533 { 534 mpq_t tmp; 535 mp_result res; 536 int out = INT_MAX; 537 538 if((res = mp_rat_init(K, &tmp)) != MP_OK) 539 return out; 540 if((res = mp_rat_set_value(K, &tmp, n, d)) != MP_OK) 541 goto CLEANUP; 542 543 out = mp_rat_compare(K, r, &tmp); 544 545 CLEANUP: 546 mp_rat_clear(K, &tmp); 547 return out; 548 } 549 550 /* }}} */ 551 552 /* {{{ mp_rat_is_integer(r) */ 553 554 int mp_rat_is_integer(mp_rat r) 555 { 556 return (mp_int_compare_value(MP_DENOM_P(r), 1) == 0); 557 } 558 559 /* }}} */ 560 561 /* {{{ mp_rat_to_ints(r, *num, *den) */ 562 563 mp_result mp_rat_to_ints(mp_rat r, mp_small *num, mp_small *den) 564 { 565 mp_result res; 566 567 if((res = mp_int_to_int(MP_NUMER_P(r), num)) != MP_OK) 568 return res; 569 570 res = mp_int_to_int(MP_DENOM_P(r), den); 571 return res; 572 } 573 574 /* }}} */ 575 576 /* {{{ mp_rat_to_string(r, radix, *str, limit) */ 577 578 mp_result mp_rat_to_string(klisp_State *K, mp_rat r, mp_size radix, char *str, 579 int limit) 580 { 581 char *start; 582 int len; 583 mp_result res; 584 585 /* Write the numerator. The sign of the rational number is written 586 by the underlying integer implementation. */ 587 if((res = mp_int_to_string(K, MP_NUMER_P(r), radix, str, limit)) != MP_OK) 588 return res; 589 590 /* If the value is zero, don't bother writing any denominator */ 591 if(mp_int_compare_zero(MP_NUMER_P(r)) == 0) 592 return MP_OK; 593 594 /* Locate the end of the numerator, and make sure we are not going to 595 exceed the limit by writing a slash. */ 596 len = strlen(str); 597 start = str + len; 598 limit -= len; 599 if(limit == 0) 600 return MP_TRUNC; 601 602 *start++ = '/'; 603 limit -= 1; 604 605 res = mp_int_to_string(K, MP_DENOM_P(r), radix, start, limit); 606 return res; 607 } 608 609 /* }}} */ 610 611 /* {{{ mp_rat_to_decimal(r, radix, prec, *str, limit) */ 612 mp_result mp_rat_to_decimal(klisp_State *K, mp_rat r, mp_size radix, 613 mp_size prec, mp_round_mode round, char *str, 614 int limit) 615 { 616 mpz_t temp[3]; 617 mp_result res; 618 char *start = str; 619 int len, lead_0, left = limit, last = 0; 620 621 SETUP(mp_int_init_copy(K, TEMP(last), MP_NUMER_P(r)), last); 622 SETUP(mp_int_init(TEMP(last)), last); 623 SETUP(mp_int_init(TEMP(last)), last); 624 625 /* Get the unsigned integer part by dividing denominator into the 626 absolute value of the numerator. */ 627 mp_int_abs(K, TEMP(0), TEMP(0)); 628 if((res = mp_int_div(K, TEMP(0), MP_DENOM_P(r), TEMP(0), 629 TEMP(1))) != MP_OK) 630 goto CLEANUP; 631 632 /* Now: T0 = integer portion, unsigned; 633 T1 = remainder, from which fractional part is computed. */ 634 635 /* Count up leading zeroes after the radix point. */ 636 for(lead_0 = 0; lead_0 < prec && mp_int_compare(TEMP(1), MP_DENOM_P(r)) < 0; 637 ++lead_0) { 638 if((res = mp_int_mul_value(K, TEMP(1), radix, TEMP(1))) != MP_OK) 639 goto CLEANUP; 640 } 641 642 /* Multiply remainder by a power of the radix sufficient to get the 643 right number of significant figures. */ 644 if(prec > lead_0) { 645 if((res = mp_int_expt_value(K, radix, prec - lead_0, 646 TEMP(2))) != MP_OK) 647 goto CLEANUP; 648 if((res = mp_int_mul(K, TEMP(1), TEMP(2), TEMP(1))) != MP_OK) 649 goto CLEANUP; 650 } 651 if((res = mp_int_div(K, TEMP(1), MP_DENOM_P(r), TEMP(1), 652 TEMP(2))) != MP_OK) 653 goto CLEANUP; 654 655 /* Now: T1 = significant digits of fractional part; 656 T2 = leftovers, to use for rounding. 657 658 At this point, what we do depends on the rounding mode. The 659 default is MP_ROUND_DOWN, for which everything is as it should be 660 already. 661 */ 662 switch(round) { 663 int cmp; 664 665 case MP_ROUND_UP: 666 if(mp_int_compare_zero(TEMP(2)) != 0) { 667 if(prec == 0) 668 res = mp_int_add_value(K, TEMP(0), 1, TEMP(0)); 669 else 670 res = mp_int_add_value(K, TEMP(1), 1, TEMP(1)); 671 } 672 break; 673 674 case MP_ROUND_HALF_UP: 675 case MP_ROUND_HALF_DOWN: 676 if((res = mp_int_mul_pow2(K, TEMP(2), 1, TEMP(2))) != MP_OK) 677 goto CLEANUP; 678 679 cmp = mp_int_compare(TEMP(2), MP_DENOM_P(r)); 680 681 if(round == MP_ROUND_HALF_UP) 682 cmp += 1; 683 684 if(cmp > 0) { 685 if(prec == 0) 686 res = mp_int_add_value(K, TEMP(0), 1, TEMP(0)); 687 else 688 res = mp_int_add_value(K, TEMP(1), 1, TEMP(1)); 689 } 690 break; 691 692 case MP_ROUND_DOWN: 693 break; /* No action required */ 694 695 default: 696 return MP_BADARG; /* Invalid rounding specifier */ 697 } 698 699 /* The sign of the output should be the sign of the numerator, but 700 if all the displayed digits will be zero due to the precision, a 701 negative shouldn't be shown. */ 702 if(MP_SIGN(MP_NUMER_P(r)) == MP_NEG && 703 (mp_int_compare_zero(TEMP(0)) != 0 || 704 mp_int_compare_zero(TEMP(1)) != 0)) { 705 *start++ = '-'; 706 left -= 1; 707 } 708 709 if((res = mp_int_to_string(K, TEMP(0), radix, start, left)) != MP_OK) 710 goto CLEANUP; 711 712 len = strlen(start); 713 start += len; 714 left -= len; 715 716 if(prec == 0) 717 goto CLEANUP; 718 719 *start++ = '.'; 720 left -= 1; 721 722 if(left < prec + 1) { 723 res = MP_TRUNC; 724 goto CLEANUP; 725 } 726 727 memset(start, '0', lead_0 - 1); 728 left -= lead_0; 729 start += lead_0 - 1; 730 731 res = mp_int_to_string(K, TEMP(1), radix, start, left); 732 733 CLEANUP: 734 while(--last >= 0) 735 mp_int_clear(K, TEMP(last)); 736 737 return res; 738 } 739 740 /* }}} */ 741 742 /* {{{ mp_rat_string_len(r, radix) */ 743 744 mp_result mp_rat_string_len(mp_rat r, mp_size radix) 745 { 746 mp_result n_len, d_len = 0; 747 748 n_len = mp_int_string_len(MP_NUMER_P(r), radix); 749 750 if(mp_int_compare_zero(MP_NUMER_P(r)) != 0) 751 d_len = mp_int_string_len(MP_DENOM_P(r), radix); 752 753 /* Though simplistic, this formula is correct. Space for the sign 754 flag is included in n_len, and the space for the NUL that is 755 counted in n_len counts for the separator here. The space for 756 the NUL counted in d_len counts for the final terminator here. */ 757 758 return n_len + d_len; 759 760 } 761 762 /* }}} */ 763 764 /* {{{ mp_rat_decimal_len(r, radix, prec) */ 765 766 mp_result mp_rat_decimal_len(mp_rat r, mp_size radix, mp_size prec) 767 { 768 int z_len, f_len; 769 770 z_len = mp_int_string_len(MP_NUMER_P(r), radix); 771 772 if(prec == 0) 773 f_len = 1; /* terminator only */ 774 else 775 f_len = 1 + prec + 1; /* decimal point, digits, terminator */ 776 777 return z_len + f_len; 778 } 779 780 /* }}} */ 781 782 /* {{{ mp_rat_read_string(r, radix, *str) */ 783 784 mp_result mp_rat_read_string(klisp_State *K, mp_rat r, mp_size radix, 785 const char *str) 786 { 787 return mp_rat_read_cstring(K, r, radix, str, NULL); 788 } 789 790 /* }}} */ 791 792 /* {{{ mp_rat_read_cstring(r, radix, *str, **end) */ 793 794 mp_result mp_rat_read_cstring(klisp_State *K, mp_rat r, mp_size radix, 795 const char *str, char **end) 796 { 797 mp_result res; 798 char *endp; 799 800 if((res = mp_int_read_cstring(K, MP_NUMER_P(r), radix, str, 801 &endp)) != MP_OK && (res != MP_TRUNC)) 802 return res; 803 804 /* Skip whitespace between numerator and (possible) separator */ 805 while(isspace((unsigned char) *endp)) 806 ++endp; 807 808 /* If there is no separator, we will stop reading at this point. */ 809 if(*endp != '/') { 810 mp_int_set_value(K, MP_DENOM_P(r), 1); 811 if(end != NULL) 812 *end = endp; 813 return res; 814 } 815 816 ++endp; /* skip separator */ 817 if((res = mp_int_read_cstring(K, MP_DENOM_P(r), radix, endp, end)) != MP_OK) 818 return res; 819 820 /* Make sure the value is well-defined */ 821 if(mp_int_compare_zero(MP_DENOM_P(r)) == 0) 822 return MP_UNDEF; 823 824 /* Reduce to lowest terms */ 825 return s_rat_reduce(K, r); 826 } 827 828 /* }}} */ 829 830 /* {{{ mp_rat_read_ustring(r, radix, *str, **end) */ 831 832 /* Read a string and figure out what format it's in. The radix may be 833 supplied as zero to use "default" behaviour. 834 835 This function will accept either a/b notation or decimal notation. 836 */ 837 mp_result mp_rat_read_ustring(klisp_State *K, mp_rat r, mp_size radix, 838 const char *str, char **end) 839 { 840 char *endp; 841 mp_result res; 842 843 if(radix == 0) 844 radix = 10; /* default to decimal input */ 845 846 if((res = mp_rat_read_cstring(K, r, radix, str, &endp)) != MP_OK) { 847 if(res == MP_TRUNC) { 848 if(*endp == '.') 849 res = mp_rat_read_cdecimal(K, r, radix, str, &endp); 850 } 851 else 852 return res; 853 } 854 855 if(end != NULL) 856 *end = endp; 857 858 return res; 859 } 860 861 /* }}} */ 862 863 /* {{{ mp_rat_read_decimal(r, radix, *str) */ 864 865 mp_result mp_rat_read_decimal(klisp_State *K, mp_rat r, mp_size radix, 866 const char *str) 867 { 868 return mp_rat_read_cdecimal(K, r, radix, str, NULL); 869 } 870 871 /* }}} */ 872 873 /* {{{ mp_rat_read_cdecimal(r, radix, *str, **end) */ 874 875 mp_result mp_rat_read_cdecimal(klisp_State *K, mp_rat r, mp_size radix, 876 const char *str, char **end) 877 { 878 mp_result res; 879 mp_sign osign; 880 char *endp; 881 882 while(isspace((unsigned char) *str)) 883 ++str; 884 885 switch(*str) { 886 case '-': 887 osign = MP_NEG; 888 break; 889 default: 890 osign = MP_ZPOS; 891 } 892 893 if((res = mp_int_read_cstring(K, MP_NUMER_P(r), radix, str, 894 &endp)) != MP_OK && (res != MP_TRUNC)) 895 return res; 896 897 /* This needs to be here. */ 898 (void) mp_int_set_value(K, MP_DENOM_P(r), 1); 899 900 if(*endp != '.') { 901 if(end != NULL) 902 *end = endp; 903 return res; 904 } 905 906 /* If the character following the decimal point is whitespace or a 907 sign flag, we will consider this a truncated value. This special 908 case is because mp_int_read_string() will consider whitespace or 909 sign flags to be valid starting characters for a value, and we do 910 not want them following the decimal point. 911 912 Once we have done this check, it is safe to read in the value of 913 the fractional piece as a regular old integer. 914 */ 915 ++endp; 916 if(*endp == '\0') { 917 if(end != NULL) 918 *end = endp; 919 return MP_OK; 920 } 921 else if(isspace((unsigned char) *endp) || *endp == '-' || *endp == '+') { 922 return MP_TRUNC; 923 } 924 else { 925 mpz_t frac; 926 mp_result save_res; 927 char *save = endp; 928 int num_lz = 0; 929 930 /* Make a temporary to hold the part after the decimal point. */ 931 if((res = mp_int_init(&frac)) != MP_OK) 932 return res; 933 934 if((res = mp_int_read_cstring(K, &frac, radix, endp, &endp)) != MP_OK && 935 (res != MP_TRUNC)) 936 goto CLEANUP; 937 938 /* Save this response for later. */ 939 save_res = res; 940 941 if(mp_int_compare_zero(&frac) == 0) 942 goto FINISHED; 943 944 /* Discard trailing zeroes (somewhat inefficiently) */ 945 while(mp_int_divisible_value(K, &frac, radix)) 946 if((res = mp_int_div_value(K, &frac, radix, &frac, NULL)) != MP_OK) 947 goto CLEANUP; 948 949 /* Count leading zeros after the decimal point */ 950 while(save[num_lz] == '0') 951 ++num_lz; 952 953 /* Find the least power of the radix that is at least as large as 954 the significant value of the fractional part, ignoring leading 955 zeroes. */ 956 (void) mp_int_set_value(K, MP_DENOM_P(r), radix); 957 958 while(mp_int_compare(MP_DENOM_P(r), &frac) < 0) { 959 if((res = mp_int_mul_value(K, MP_DENOM_P(r), radix, 960 MP_DENOM_P(r))) != MP_OK) 961 goto CLEANUP; 962 } 963 964 /* Also shift by enough to account for leading zeroes */ 965 while(num_lz > 0) { 966 if((res = mp_int_mul_value(K, MP_DENOM_P(r), radix, 967 MP_DENOM_P(r))) != MP_OK) 968 goto CLEANUP; 969 970 --num_lz; 971 } 972 973 /* Having found this power, shift the numerator leftward that 974 many, digits, and add the nonzero significant digits of the 975 fractional part to get the result. */ 976 if((res = mp_int_mul(K, MP_NUMER_P(r), MP_DENOM_P(r), 977 MP_NUMER_P(r))) != MP_OK) 978 goto CLEANUP; 979 980 { /* This addition needs to be unsigned. */ 981 MP_SIGN(MP_NUMER_P(r)) = MP_ZPOS; 982 if((res = mp_int_add(K, MP_NUMER_P(r), &frac, 983 MP_NUMER_P(r))) != MP_OK) 984 goto CLEANUP; 985 986 MP_SIGN(MP_NUMER_P(r)) = osign; 987 } 988 if((res = s_rat_reduce(K, r)) != MP_OK) 989 goto CLEANUP; 990 991 /* At this point, what we return depends on whether reading the 992 fractional part was truncated or not. That information is 993 saved from when we called mp_int_read_string() above. */ 994 FINISHED: 995 res = save_res; 996 if(end != NULL) 997 *end = endp; 998 999 CLEANUP: 1000 mp_int_clear(K, &frac); 1001 1002 return res; 1003 } 1004 } 1005 1006 /* }}} */ 1007 1008 /* Private functions for internal use. Make unchecked assumptions 1009 about format and validity of inputs. */ 1010 1011 /* {{{ s_rat_reduce(r) */ 1012 1013 static mp_result s_rat_reduce(klisp_State *K, mp_rat r) 1014 { 1015 mpz_t gcd; 1016 mp_result res = MP_OK; 1017 1018 if(mp_int_compare_zero(MP_NUMER_P(r)) == 0) { 1019 mp_int_set_value(K, MP_DENOM_P(r), 1); 1020 return MP_OK; 1021 } 1022 1023 /* If the greatest common divisor of the numerator and denominator 1024 is greater than 1, divide it out. */ 1025 if((res = mp_int_init(&gcd)) != MP_OK) 1026 return res; 1027 1028 if((res = mp_int_gcd(K, MP_NUMER_P(r), MP_DENOM_P(r), &gcd)) != MP_OK) 1029 goto CLEANUP; 1030 1031 if(mp_int_compare_value(&gcd, 1) != 0) { 1032 if((res = mp_int_div(K, MP_NUMER_P(r), &gcd, MP_NUMER_P(r), 1033 NULL)) != MP_OK) 1034 goto CLEANUP; 1035 if((res = mp_int_div(K, MP_DENOM_P(r), &gcd, MP_DENOM_P(r), 1036 NULL)) != MP_OK) 1037 goto CLEANUP; 1038 } 1039 1040 /* Fix up the signs of numerator and denominator */ 1041 if(MP_SIGN(MP_NUMER_P(r)) == MP_SIGN(MP_DENOM_P(r))) 1042 MP_SIGN(MP_NUMER_P(r)) = MP_SIGN(MP_DENOM_P(r)) = MP_ZPOS; 1043 else { 1044 MP_SIGN(MP_NUMER_P(r)) = MP_NEG; 1045 MP_SIGN(MP_DENOM_P(r)) = MP_ZPOS; 1046 } 1047 1048 CLEANUP: 1049 mp_int_clear(K, &gcd); 1050 1051 return res; 1052 } 1053 1054 /* }}} */ 1055 1056 /* {{{ s_rat_combine(a, b, c, comb_f) */ 1057 1058 static mp_result s_rat_combine(klisp_State *K, mp_rat a, mp_rat b, mp_rat c, 1059 mp_result (*comb_f)(klisp_State *K, mp_int, 1060 mp_int, mp_int)) 1061 { 1062 mp_result res; 1063 1064 /* Shortcut when denominators are already common */ 1065 if(mp_int_compare(MP_DENOM_P(a), MP_DENOM_P(b)) == 0) { 1066 if((res = (comb_f)(K, MP_NUMER_P(a), MP_NUMER_P(b), 1067 MP_NUMER_P(c))) != MP_OK) 1068 return res; 1069 if((res = mp_int_copy(K, MP_DENOM_P(a), MP_DENOM_P(c))) != MP_OK) 1070 return res; 1071 1072 return s_rat_reduce(K, c); 1073 } 1074 else { 1075 mpz_t temp[2]; 1076 int last = 0; 1077 1078 SETUP(mp_int_init_copy(K, TEMP(last), MP_NUMER_P(a)), last); 1079 SETUP(mp_int_init_copy(K, TEMP(last), MP_NUMER_P(b)), last); 1080 1081 if((res = mp_int_mul(K, TEMP(0), MP_DENOM_P(b), TEMP(0))) != MP_OK) 1082 goto CLEANUP; 1083 if((res = mp_int_mul(K, TEMP(1), MP_DENOM_P(a), TEMP(1))) != MP_OK) 1084 goto CLEANUP; 1085 if((res = (comb_f)(K, TEMP(0), TEMP(1), MP_NUMER_P(c))) != MP_OK) 1086 goto CLEANUP; 1087 1088 res = mp_int_mul(K, MP_DENOM_P(a), MP_DENOM_P(b), MP_DENOM_P(c)); 1089 1090 CLEANUP: 1091 while(--last >= 0) 1092 mp_int_clear(K, TEMP(last)); 1093 1094 if(res == MP_OK) 1095 return s_rat_reduce(K, c); 1096 else 1097 return res; 1098 } 1099 } 1100 1101 /* }}} */ 1102 1103 /* Here there be dragons */