941 lines · c
1/*2 Name: imrat.c3 Purpose: Arbitrary precision rational arithmetic routines.4 Author: M. J. Fromberger5 6 Copyright (C) 2002-2007 Michael J. Fromberger, All Rights Reserved.7 8 Permission is hereby granted, free of charge, to any person obtaining a copy9 of this software and associated documentation files (the "Software"), to deal10 in the Software without restriction, including without limitation the rights11 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell12 copies of the Software, and to permit persons to whom the Software is13 furnished to do so, subject to the following conditions:14 15 The above copyright notice and this permission notice shall be included in16 all copies or substantial portions of the Software.17 18 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR19 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,20 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE21 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER22 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,23 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE24 SOFTWARE.25 */26 27#include "imrat.h"28#include <assert.h>29#include <ctype.h>30#include <stdlib.h>31#include <string.h>32 33#define MP_NUMER_SIGN(Q) (MP_NUMER_P(Q)->sign)34#define MP_DENOM_SIGN(Q) (MP_DENOM_P(Q)->sign)35 36#define TEMP(K) (temp + (K))37#define SETUP(E, C) \38 do { \39 if ((res = (E)) != MP_OK) goto CLEANUP; \40 ++(C); \41 } while (0)42 43/* Reduce the given rational, in place, to lowest terms and canonical form.44 Zero is represented as 0/1, one as 1/1. Signs are adjusted so that the sign45 of the numerator is definitive. */46static mp_result s_rat_reduce(mp_rat r);47 48/* Common code for addition and subtraction operations on rationals. */49static mp_result s_rat_combine(mp_rat a, mp_rat b, mp_rat c,50 mp_result (*comb_f)(mp_int, mp_int, mp_int));51 52mp_result mp_rat_init(mp_rat r) {53 mp_result res;54 55 if ((res = mp_int_init(MP_NUMER_P(r))) != MP_OK) return res;56 if ((res = mp_int_init(MP_DENOM_P(r))) != MP_OK) {57 mp_int_clear(MP_NUMER_P(r));58 return res;59 }60 61 return mp_int_set_value(MP_DENOM_P(r), 1);62}63 64mp_rat mp_rat_alloc(void) {65 mp_rat out = malloc(sizeof(*out));66 67 if (out != NULL) {68 if (mp_rat_init(out) != MP_OK) {69 free(out);70 return NULL;71 }72 }73 74 return out;75}76 77mp_result mp_rat_reduce(mp_rat r) { return s_rat_reduce(r); }78 79mp_result mp_rat_init_size(mp_rat r, mp_size n_prec, mp_size d_prec) {80 mp_result res;81 82 if ((res = mp_int_init_size(MP_NUMER_P(r), n_prec)) != MP_OK) {83 return res;84 }85 if ((res = mp_int_init_size(MP_DENOM_P(r), d_prec)) != MP_OK) {86 mp_int_clear(MP_NUMER_P(r));87 return res;88 }89 90 return mp_int_set_value(MP_DENOM_P(r), 1);91}92 93mp_result mp_rat_init_copy(mp_rat r, mp_rat old) {94 mp_result res;95 96 if ((res = mp_int_init_copy(MP_NUMER_P(r), MP_NUMER_P(old))) != MP_OK) {97 return res;98 }99 if ((res = mp_int_init_copy(MP_DENOM_P(r), MP_DENOM_P(old))) != MP_OK)100 mp_int_clear(MP_NUMER_P(r));101 102 return res;103}104 105mp_result mp_rat_set_value(mp_rat r, mp_small numer, mp_small denom) {106 mp_result res;107 108 if (denom == 0) return MP_UNDEF;109 110 if ((res = mp_int_set_value(MP_NUMER_P(r), numer)) != MP_OK) {111 return res;112 }113 if ((res = mp_int_set_value(MP_DENOM_P(r), denom)) != MP_OK) {114 return res;115 }116 117 return s_rat_reduce(r);118}119 120mp_result mp_rat_set_uvalue(mp_rat r, mp_usmall numer, mp_usmall denom) {121 mp_result res;122 123 if (denom == 0) return MP_UNDEF;124 125 if ((res = mp_int_set_uvalue(MP_NUMER_P(r), numer)) != MP_OK) {126 return res;127 }128 if ((res = mp_int_set_uvalue(MP_DENOM_P(r), denom)) != MP_OK) {129 return res;130 }131 132 return s_rat_reduce(r);133}134 135void mp_rat_clear(mp_rat r) {136 mp_int_clear(MP_NUMER_P(r));137 mp_int_clear(MP_DENOM_P(r));138}139 140void mp_rat_free(mp_rat r) {141 assert(r != NULL);142 143 if (r->num.digits != NULL) mp_rat_clear(r);144 145 free(r);146}147 148mp_result mp_rat_numer(mp_rat r, mp_int z) {149 return mp_int_copy(MP_NUMER_P(r), z);150}151 152mp_int mp_rat_numer_ref(mp_rat r) { return MP_NUMER_P(r); }153 154mp_result mp_rat_denom(mp_rat r, mp_int z) {155 return mp_int_copy(MP_DENOM_P(r), z);156}157 158mp_int mp_rat_denom_ref(mp_rat r) { return MP_DENOM_P(r); }159 160mp_sign mp_rat_sign(mp_rat r) { return MP_NUMER_SIGN(r); }161 162mp_result mp_rat_copy(mp_rat a, mp_rat c) {163 mp_result res;164 165 if ((res = mp_int_copy(MP_NUMER_P(a), MP_NUMER_P(c))) != MP_OK) {166 return res;167 }168 169 res = mp_int_copy(MP_DENOM_P(a), MP_DENOM_P(c));170 return res;171}172 173void mp_rat_zero(mp_rat r) {174 mp_int_zero(MP_NUMER_P(r));175 mp_int_set_value(MP_DENOM_P(r), 1);176}177 178mp_result mp_rat_abs(mp_rat a, mp_rat c) {179 mp_result res;180 181 if ((res = mp_int_abs(MP_NUMER_P(a), MP_NUMER_P(c))) != MP_OK) {182 return res;183 }184 185 res = mp_int_abs(MP_DENOM_P(a), MP_DENOM_P(c));186 return res;187}188 189mp_result mp_rat_neg(mp_rat a, mp_rat c) {190 mp_result res;191 192 if ((res = mp_int_neg(MP_NUMER_P(a), MP_NUMER_P(c))) != MP_OK) {193 return res;194 }195 196 res = mp_int_copy(MP_DENOM_P(a), MP_DENOM_P(c));197 return res;198}199 200mp_result mp_rat_recip(mp_rat a, mp_rat c) {201 mp_result res;202 203 if (mp_rat_compare_zero(a) == 0) return MP_UNDEF;204 205 if ((res = mp_rat_copy(a, c)) != MP_OK) return res;206 207 mp_int_swap(MP_NUMER_P(c), MP_DENOM_P(c));208 209 /* Restore the signs of the swapped elements */210 {211 mp_sign tmp = MP_NUMER_SIGN(c);212 213 MP_NUMER_SIGN(c) = MP_DENOM_SIGN(c);214 MP_DENOM_SIGN(c) = tmp;215 }216 217 return MP_OK;218}219 220mp_result mp_rat_add(mp_rat a, mp_rat b, mp_rat c) {221 return s_rat_combine(a, b, c, mp_int_add);222}223 224mp_result mp_rat_sub(mp_rat a, mp_rat b, mp_rat c) {225 return s_rat_combine(a, b, c, mp_int_sub);226}227 228mp_result mp_rat_mul(mp_rat a, mp_rat b, mp_rat c) {229 mp_result res;230 231 if ((res = mp_int_mul(MP_NUMER_P(a), MP_NUMER_P(b), MP_NUMER_P(c))) != MP_OK)232 return res;233 234 if (mp_int_compare_zero(MP_NUMER_P(c)) != 0) {235 res = mp_int_mul(MP_DENOM_P(a), MP_DENOM_P(b), MP_DENOM_P(c));236 if (res != MP_OK) {237 return res;238 }239 }240 241 return s_rat_reduce(c);242}243 244mp_result mp_rat_div(mp_rat a, mp_rat b, mp_rat c) {245 mp_result res = MP_OK;246 247 if (mp_rat_compare_zero(b) == 0) return MP_UNDEF;248 249 if (c == a || c == b) {250 mpz_t tmp;251 252 if ((res = mp_int_init(&tmp)) != MP_OK) return res;253 if ((res = mp_int_mul(MP_NUMER_P(a), MP_DENOM_P(b), &tmp)) != MP_OK) {254 goto CLEANUP;255 }256 if ((res = mp_int_mul(MP_DENOM_P(a), MP_NUMER_P(b), MP_DENOM_P(c))) !=257 MP_OK) {258 goto CLEANUP;259 }260 res = mp_int_copy(&tmp, MP_NUMER_P(c));261 262 CLEANUP:263 mp_int_clear(&tmp);264 } else {265 if ((res = mp_int_mul(MP_NUMER_P(a), MP_DENOM_P(b), MP_NUMER_P(c))) !=266 MP_OK) {267 return res;268 }269 if ((res = mp_int_mul(MP_DENOM_P(a), MP_NUMER_P(b), MP_DENOM_P(c))) !=270 MP_OK) {271 return res;272 }273 }274 275 if (res != MP_OK) {276 return res;277 } else {278 return s_rat_reduce(c);279 }280}281 282mp_result mp_rat_add_int(mp_rat a, mp_int b, mp_rat c) {283 mpz_t tmp;284 mp_result res;285 286 if ((res = mp_int_init_copy(&tmp, b)) != MP_OK) {287 return res;288 }289 if ((res = mp_int_mul(&tmp, MP_DENOM_P(a), &tmp)) != MP_OK) {290 goto CLEANUP;291 }292 if ((res = mp_rat_copy(a, c)) != MP_OK) {293 goto CLEANUP;294 }295 if ((res = mp_int_add(MP_NUMER_P(c), &tmp, MP_NUMER_P(c))) != MP_OK) {296 goto CLEANUP;297 }298 299 res = s_rat_reduce(c);300 301CLEANUP:302 mp_int_clear(&tmp);303 return res;304}305 306mp_result mp_rat_sub_int(mp_rat a, mp_int b, mp_rat c) {307 mpz_t tmp;308 mp_result res;309 310 if ((res = mp_int_init_copy(&tmp, b)) != MP_OK) {311 return res;312 }313 if ((res = mp_int_mul(&tmp, MP_DENOM_P(a), &tmp)) != MP_OK) {314 goto CLEANUP;315 }316 if ((res = mp_rat_copy(a, c)) != MP_OK) {317 goto CLEANUP;318 }319 if ((res = mp_int_sub(MP_NUMER_P(c), &tmp, MP_NUMER_P(c))) != MP_OK) {320 goto CLEANUP;321 }322 323 res = s_rat_reduce(c);324 325CLEANUP:326 mp_int_clear(&tmp);327 return res;328}329 330mp_result mp_rat_mul_int(mp_rat a, mp_int b, mp_rat c) {331 mp_result res;332 333 if ((res = mp_rat_copy(a, c)) != MP_OK) {334 return res;335 }336 if ((res = mp_int_mul(MP_NUMER_P(c), b, MP_NUMER_P(c))) != MP_OK) {337 return res;338 }339 340 return s_rat_reduce(c);341}342 343mp_result mp_rat_div_int(mp_rat a, mp_int b, mp_rat c) {344 mp_result res;345 346 if (mp_int_compare_zero(b) == 0) {347 return MP_UNDEF;348 }349 if ((res = mp_rat_copy(a, c)) != MP_OK) {350 return res;351 }352 if ((res = mp_int_mul(MP_DENOM_P(c), b, MP_DENOM_P(c))) != MP_OK) {353 return res;354 }355 356 return s_rat_reduce(c);357}358 359mp_result mp_rat_expt(mp_rat a, mp_small b, mp_rat c) {360 mp_result res;361 362 /* Special cases for easy powers. */363 if (b == 0) {364 return mp_rat_set_value(c, 1, 1);365 } else if (b == 1) {366 return mp_rat_copy(a, c);367 }368 369 /* Since rationals are always stored in lowest terms, it is not necessary to370 reduce again when raising to an integer power. */371 if ((res = mp_int_expt(MP_NUMER_P(a), b, MP_NUMER_P(c))) != MP_OK) {372 return res;373 }374 375 return mp_int_expt(MP_DENOM_P(a), b, MP_DENOM_P(c));376}377 378int mp_rat_compare(mp_rat a, mp_rat b) {379 /* Quick check for opposite signs. Works because the sign of the numerator380 is always definitive. */381 if (MP_NUMER_SIGN(a) != MP_NUMER_SIGN(b)) {382 if (MP_NUMER_SIGN(a) == MP_ZPOS) {383 return 1;384 } else {385 return -1;386 }387 } else {388 /* Compare absolute magnitudes; if both are positive, the answer stands,389 otherwise it needs to be reflected about zero. */390 int cmp = mp_rat_compare_unsigned(a, b);391 392 if (MP_NUMER_SIGN(a) == MP_ZPOS) {393 return cmp;394 } else {395 return -cmp;396 }397 }398}399 400int mp_rat_compare_unsigned(mp_rat a, mp_rat b) {401 /* If the denominators are equal, we can quickly compare numerators without402 multiplying. Otherwise, we actually have to do some work. */403 if (mp_int_compare_unsigned(MP_DENOM_P(a), MP_DENOM_P(b)) == 0) {404 return mp_int_compare_unsigned(MP_NUMER_P(a), MP_NUMER_P(b));405 }406 407 else {408 mpz_t temp[2];409 mp_result res;410 int cmp = INT_MAX, last = 0;411 412 /* t0 = num(a) * den(b), t1 = num(b) * den(a) */413 SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(a)), last);414 SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(b)), last);415 416 if ((res = mp_int_mul(TEMP(0), MP_DENOM_P(b), TEMP(0))) != MP_OK ||417 (res = mp_int_mul(TEMP(1), MP_DENOM_P(a), TEMP(1))) != MP_OK) {418 goto CLEANUP;419 }420 421 cmp = mp_int_compare_unsigned(TEMP(0), TEMP(1));422 423 CLEANUP:424 while (--last >= 0) {425 mp_int_clear(TEMP(last));426 }427 428 return cmp;429 }430}431 432int mp_rat_compare_zero(mp_rat r) { return mp_int_compare_zero(MP_NUMER_P(r)); }433 434int mp_rat_compare_value(mp_rat r, mp_small n, mp_small d) {435 mpq_t tmp;436 mp_result res;437 int out = INT_MAX;438 439 if ((res = mp_rat_init(&tmp)) != MP_OK) {440 return out;441 }442 if ((res = mp_rat_set_value(&tmp, n, d)) != MP_OK) {443 goto CLEANUP;444 }445 446 out = mp_rat_compare(r, &tmp);447 448CLEANUP:449 mp_rat_clear(&tmp);450 return out;451}452 453bool mp_rat_is_integer(mp_rat r) {454 return (mp_int_compare_value(MP_DENOM_P(r), 1) == 0);455}456 457mp_result mp_rat_to_ints(mp_rat r, mp_small *num, mp_small *den) {458 mp_result res;459 460 if ((res = mp_int_to_int(MP_NUMER_P(r), num)) != MP_OK) {461 return res;462 }463 464 res = mp_int_to_int(MP_DENOM_P(r), den);465 return res;466}467 468mp_result mp_rat_to_string(mp_rat r, mp_size radix, char *str, int limit) {469 /* Write the numerator. The sign of the rational number is written by the470 underlying integer implementation. */471 mp_result res;472 if ((res = mp_int_to_string(MP_NUMER_P(r), radix, str, limit)) != MP_OK) {473 return res;474 }475 476 /* If the value is zero, don't bother writing any denominator */477 if (mp_int_compare_zero(MP_NUMER_P(r)) == 0) {478 return MP_OK;479 }480 481 /* Locate the end of the numerator, and make sure we are not going to exceed482 the limit by writing a slash. */483 int len = strlen(str);484 char *start = str + len;485 limit -= len;486 if (limit == 0) return MP_TRUNC;487 488 *start++ = '/';489 limit -= 1;490 491 return mp_int_to_string(MP_DENOM_P(r), radix, start, limit);492}493 494mp_result mp_rat_to_decimal(mp_rat r, mp_size radix, mp_size prec,495 mp_round_mode round, char *str, int limit) {496 mpz_t temp[3];497 mp_result res;498 int last = 0;499 500 SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(r)), last);501 SETUP(mp_int_init(TEMP(last)), last);502 SETUP(mp_int_init(TEMP(last)), last);503 504 /* Get the unsigned integer part by dividing denominator into the absolute505 value of the numerator. */506 mp_int_abs(TEMP(0), TEMP(0));507 if ((res = mp_int_div(TEMP(0), MP_DENOM_P(r), TEMP(0), TEMP(1))) != MP_OK) {508 goto CLEANUP;509 }510 511 /* Now: T0 = integer portion, unsigned;512 T1 = remainder, from which fractional part is computed. */513 514 /* Count up leading zeroes after the radix point. */515 int zprec = (int)prec;516 int lead_0;517 for (lead_0 = 0; lead_0 < zprec && mp_int_compare(TEMP(1), MP_DENOM_P(r)) < 0;518 ++lead_0) {519 if ((res = mp_int_mul_value(TEMP(1), radix, TEMP(1))) != MP_OK) {520 goto CLEANUP;521 }522 }523 524 /* Multiply remainder by a power of the radix sufficient to get the right525 number of significant figures. */526 if (zprec > lead_0) {527 if ((res = mp_int_expt_value(radix, zprec - lead_0, TEMP(2))) != MP_OK) {528 goto CLEANUP;529 }530 if ((res = mp_int_mul(TEMP(1), TEMP(2), TEMP(1))) != MP_OK) {531 goto CLEANUP;532 }533 }534 if ((res = mp_int_div(TEMP(1), MP_DENOM_P(r), TEMP(1), TEMP(2))) != MP_OK) {535 goto CLEANUP;536 }537 538 /* Now: T1 = significant digits of fractional part;539 T2 = leftovers, to use for rounding.540 541 At this point, what we do depends on the rounding mode. The default is542 MP_ROUND_DOWN, for which everything is as it should be already.543 */544 switch (round) {545 int cmp;546 547 case MP_ROUND_UP:548 if (mp_int_compare_zero(TEMP(2)) != 0) {549 if (prec == 0) {550 res = mp_int_add_value(TEMP(0), 1, TEMP(0));551 } else {552 res = mp_int_add_value(TEMP(1), 1, TEMP(1));553 }554 }555 break;556 557 case MP_ROUND_HALF_UP:558 case MP_ROUND_HALF_DOWN:559 if ((res = mp_int_mul_pow2(TEMP(2), 1, TEMP(2))) != MP_OK) {560 goto CLEANUP;561 }562 563 cmp = mp_int_compare(TEMP(2), MP_DENOM_P(r));564 565 if (round == MP_ROUND_HALF_UP) cmp += 1;566 567 if (cmp > 0) {568 if (prec == 0) {569 res = mp_int_add_value(TEMP(0), 1, TEMP(0));570 } else {571 res = mp_int_add_value(TEMP(1), 1, TEMP(1));572 }573 }574 break;575 576 case MP_ROUND_DOWN:577 break; /* No action required */578 579 default:580 return MP_BADARG; /* Invalid rounding specifier */581 }582 if (res != MP_OK) {583 goto CLEANUP;584 }585 586 /* The sign of the output should be the sign of the numerator, but if all the587 displayed digits will be zero due to the precision, a negative shouldn't588 be shown. */589 char *start = str;590 int left = limit;591 if (MP_NUMER_SIGN(r) == MP_NEG && (mp_int_compare_zero(TEMP(0)) != 0 ||592 mp_int_compare_zero(TEMP(1)) != 0)) {593 *start++ = '-';594 left -= 1;595 }596 597 if ((res = mp_int_to_string(TEMP(0), radix, start, left)) != MP_OK) {598 goto CLEANUP;599 }600 601 int len = strlen(start);602 start += len;603 left -= len;604 605 if (prec == 0) goto CLEANUP;606 607 *start++ = '.';608 left -= 1;609 610 if (left < zprec + 1) {611 res = MP_TRUNC;612 goto CLEANUP;613 }614 615 memset(start, '0', lead_0 - 1);616 left -= lead_0;617 start += lead_0 - 1;618 619 res = mp_int_to_string(TEMP(1), radix, start, left);620 621CLEANUP:622 while (--last >= 0) mp_int_clear(TEMP(last));623 624 return res;625}626 627mp_result mp_rat_string_len(mp_rat r, mp_size radix) {628 mp_result d_len = 0;629 mp_result n_len = mp_int_string_len(MP_NUMER_P(r), radix);630 631 if (mp_int_compare_zero(MP_NUMER_P(r)) != 0) {632 d_len = mp_int_string_len(MP_DENOM_P(r), radix);633 }634 635 /* Though simplistic, this formula is correct. Space for the sign flag is636 included in n_len, and the space for the NUL that is counted in n_len637 counts for the separator here. The space for the NUL counted in d_len638 counts for the final terminator here. */639 640 return n_len + d_len;641}642 643mp_result mp_rat_decimal_len(mp_rat r, mp_size radix, mp_size prec) {644 int f_len;645 int z_len = mp_int_string_len(MP_NUMER_P(r), radix);646 647 if (prec == 0) {648 f_len = 1; /* terminator only */649 } else {650 f_len = 1 + prec + 1; /* decimal point, digits, terminator */651 }652 653 return z_len + f_len;654}655 656mp_result mp_rat_read_string(mp_rat r, mp_size radix, const char *str) {657 return mp_rat_read_cstring(r, radix, str, NULL);658}659 660mp_result mp_rat_read_cstring(mp_rat r, mp_size radix, const char *str,661 char **end) {662 mp_result res;663 char *endp;664 665 if ((res = mp_int_read_cstring(MP_NUMER_P(r), radix, str, &endp)) != MP_OK &&666 (res != MP_TRUNC)) {667 return res;668 }669 670 /* Skip whitespace between numerator and (possible) separator */671 while (isspace((unsigned char)*endp)) {672 ++endp;673 }674 675 /* If there is no separator, we will stop reading at this point. */676 if (*endp != '/') {677 mp_int_set_value(MP_DENOM_P(r), 1);678 if (end != NULL) *end = endp;679 return res;680 }681 682 ++endp; /* skip separator */683 if ((res = mp_int_read_cstring(MP_DENOM_P(r), radix, endp, end)) != MP_OK) {684 return res;685 }686 687 /* Make sure the value is well-defined */688 if (mp_int_compare_zero(MP_DENOM_P(r)) == 0) {689 return MP_UNDEF;690 }691 692 /* Reduce to lowest terms */693 return s_rat_reduce(r);694}695 696/* Read a string and figure out what format it's in. The radix may be supplied697 as zero to use "default" behaviour.698 699 This function will accept either a/b notation or decimal notation.700 */701mp_result mp_rat_read_ustring(mp_rat r, mp_size radix, const char *str,702 char **end) {703 char *endp = "";704 mp_result res;705 706 if (radix == 0) radix = 10; /* default to decimal input */707 708 res = mp_rat_read_cstring(r, radix, str, &endp);709 if (res == MP_TRUNC && *endp == '.') {710 res = mp_rat_read_cdecimal(r, radix, str, &endp);711 }712 if (end != NULL) *end = endp;713 return res;714}715 716mp_result mp_rat_read_decimal(mp_rat r, mp_size radix, const char *str) {717 return mp_rat_read_cdecimal(r, radix, str, NULL);718}719 720mp_result mp_rat_read_cdecimal(mp_rat r, mp_size radix, const char *str,721 char **end) {722 mp_result res;723 mp_sign osign;724 char *endp;725 726 while (isspace((unsigned char)*str)) ++str;727 728 switch (*str) {729 case '-':730 osign = MP_NEG;731 break;732 default:733 osign = MP_ZPOS;734 }735 736 if ((res = mp_int_read_cstring(MP_NUMER_P(r), radix, str, &endp)) != MP_OK &&737 (res != MP_TRUNC)) {738 return res;739 }740 741 /* This needs to be here. */742 (void)mp_int_set_value(MP_DENOM_P(r), 1);743 744 if (*endp != '.') {745 if (end != NULL) *end = endp;746 return res;747 }748 749 /* If the character following the decimal point is whitespace or a sign flag,750 we will consider this a truncated value. This special case is because751 mp_int_read_string() will consider whitespace or sign flags to be valid752 starting characters for a value, and we do not want them following the753 decimal point.754 755 Once we have done this check, it is safe to read in the value of the756 fractional piece as a regular old integer.757 */758 ++endp;759 if (*endp == '\0') {760 if (end != NULL) *end = endp;761 return MP_OK;762 } else if (isspace((unsigned char)*endp) || *endp == '-' || *endp == '+') {763 return MP_TRUNC;764 } else {765 mpz_t frac;766 mp_result save_res;767 char *save = endp;768 int num_lz = 0;769 770 /* Make a temporary to hold the part after the decimal point. */771 if ((res = mp_int_init(&frac)) != MP_OK) {772 return res;773 }774 775 if ((res = mp_int_read_cstring(&frac, radix, endp, &endp)) != MP_OK &&776 (res != MP_TRUNC)) {777 goto CLEANUP;778 }779 780 /* Save this response for later. */781 save_res = res;782 783 if (mp_int_compare_zero(&frac) == 0) goto FINISHED;784 785 /* Discard trailing zeroes (somewhat inefficiently) */786 while (mp_int_divisible_value(&frac, radix)) {787 if ((res = mp_int_div_value(&frac, radix, &frac, NULL)) != MP_OK) {788 goto CLEANUP;789 }790 }791 792 /* Count leading zeros after the decimal point */793 while (save[num_lz] == '0') {794 ++num_lz;795 }796 797 /* Find the least power of the radix that is at least as large as the798 significant value of the fractional part, ignoring leading zeroes. */799 (void)mp_int_set_value(MP_DENOM_P(r), radix);800 801 while (mp_int_compare(MP_DENOM_P(r), &frac) < 0) {802 if ((res = mp_int_mul_value(MP_DENOM_P(r), radix, MP_DENOM_P(r))) !=803 MP_OK) {804 goto CLEANUP;805 }806 }807 808 /* Also shift by enough to account for leading zeroes */809 while (num_lz > 0) {810 if ((res = mp_int_mul_value(MP_DENOM_P(r), radix, MP_DENOM_P(r))) !=811 MP_OK) {812 goto CLEANUP;813 }814 815 --num_lz;816 }817 818 /* Having found this power, shift the numerator leftward that many, digits,819 and add the nonzero significant digits of the fractional part to get the820 result. */821 if ((res = mp_int_mul(MP_NUMER_P(r), MP_DENOM_P(r), MP_NUMER_P(r))) !=822 MP_OK) {823 goto CLEANUP;824 }825 826 { /* This addition needs to be unsigned. */827 MP_NUMER_SIGN(r) = MP_ZPOS;828 if ((res = mp_int_add(MP_NUMER_P(r), &frac, MP_NUMER_P(r))) != MP_OK) {829 goto CLEANUP;830 }831 832 MP_NUMER_SIGN(r) = osign;833 }834 if ((res = s_rat_reduce(r)) != MP_OK) goto CLEANUP;835 836 /* At this point, what we return depends on whether reading the fractional837 part was truncated or not. That information is saved from when we838 called mp_int_read_string() above. */839 FINISHED:840 res = save_res;841 if (end != NULL) *end = endp;842 843 CLEANUP:844 mp_int_clear(&frac);845 846 return res;847 }848}849 850/* Private functions for internal use. Make unchecked assumptions about format851 and validity of inputs. */852 853static mp_result s_rat_reduce(mp_rat r) {854 mpz_t gcd;855 mp_result res = MP_OK;856 857 if (mp_int_compare_zero(MP_NUMER_P(r)) == 0) {858 mp_int_set_value(MP_DENOM_P(r), 1);859 return MP_OK;860 }861 862 /* If the greatest common divisor of the numerator and denominator is greater863 than 1, divide it out. */864 if ((res = mp_int_init(&gcd)) != MP_OK) return res;865 866 if ((res = mp_int_gcd(MP_NUMER_P(r), MP_DENOM_P(r), &gcd)) != MP_OK) {867 goto CLEANUP;868 }869 870 if (mp_int_compare_value(&gcd, 1) != 0) {871 if ((res = mp_int_div(MP_NUMER_P(r), &gcd, MP_NUMER_P(r), NULL)) != MP_OK) {872 goto CLEANUP;873 }874 if ((res = mp_int_div(MP_DENOM_P(r), &gcd, MP_DENOM_P(r), NULL)) != MP_OK) {875 goto CLEANUP;876 }877 }878 879 /* Fix up the signs of numerator and denominator */880 if (MP_NUMER_SIGN(r) == MP_DENOM_SIGN(r)) {881 MP_NUMER_SIGN(r) = MP_DENOM_SIGN(r) = MP_ZPOS;882 } else {883 MP_NUMER_SIGN(r) = MP_NEG;884 MP_DENOM_SIGN(r) = MP_ZPOS;885 }886 887CLEANUP:888 mp_int_clear(&gcd);889 890 return res;891}892 893static mp_result s_rat_combine(mp_rat a, mp_rat b, mp_rat c,894 mp_result (*comb_f)(mp_int, mp_int, mp_int)) {895 mp_result res;896 897 /* Shortcut when denominators are already common */898 if (mp_int_compare(MP_DENOM_P(a), MP_DENOM_P(b)) == 0) {899 if ((res = (comb_f)(MP_NUMER_P(a), MP_NUMER_P(b), MP_NUMER_P(c))) !=900 MP_OK) {901 return res;902 }903 if ((res = mp_int_copy(MP_DENOM_P(a), MP_DENOM_P(c))) != MP_OK) {904 return res;905 }906 907 return s_rat_reduce(c);908 } else {909 mpz_t temp[2];910 int last = 0;911 912 SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(a)), last);913 SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(b)), last);914 915 if ((res = mp_int_mul(TEMP(0), MP_DENOM_P(b), TEMP(0))) != MP_OK) {916 goto CLEANUP;917 }918 if ((res = mp_int_mul(TEMP(1), MP_DENOM_P(a), TEMP(1))) != MP_OK) {919 goto CLEANUP;920 }921 if ((res = (comb_f)(TEMP(0), TEMP(1), MP_NUMER_P(c))) != MP_OK) {922 goto CLEANUP;923 }924 925 res = mp_int_mul(MP_DENOM_P(a), MP_DENOM_P(b), MP_DENOM_P(c));926 927 CLEANUP:928 while (--last >= 0) {929 mp_int_clear(TEMP(last));930 }931 932 if (res == MP_OK) {933 return s_rat_reduce(c);934 } else {935 return res;936 }937 }938}939 940/* Here there be dragons */941