polly/lib/External/isl/imath/imrat.c

   1 /*
   2   Name:     imrat.c
   3   Purpose:  Arbitrary precision rational arithmetic routines.
   4   Author:   M. J. Fromberger
   5
   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 copy
   9   of this software and associated documentation files (the "Software"), to deal
  10   in the Software without restriction, including without limitation the rights
  11   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  12   copies of the Software, and to permit persons to whom the Software is
  13   furnished to do so, subject to the following conditions:
  14
  15   The above copyright notice and this permission notice shall be included in
  16   all copies or substantial portions of the Software.
  17
  18   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  19   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  20   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
  21   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  22   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 THE
  24   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 sign
  45    of the numerator is definitive. */
  46 static mp_result s_rat_reduce(mp_rat r);
  47
  48 /* Common code for addition and subtraction operations on rationals. */
  49 static 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
  52 mp_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
  64 mp_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
  77 mp_result mp_rat_reduce(mp_rat r) { return s_rat_reduce(r); }
  78
  79 mp_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
  93 mp_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
 105 mp_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
 120 mp_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
 135 void mp_rat_clear(mp_rat r) {
 136   mp_int_clear(MP_NUMER_P(r));
 137   mp_int_clear(MP_DENOM_P(r));
 138 }
 139
 140 void 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
 148 mp_result mp_rat_numer(mp_rat r, mp_int z) {
 149   return mp_int_copy(MP_NUMER_P(r), z);
 150 }
 151
 152 mp_int mp_rat_numer_ref(mp_rat r) { return MP_NUMER_P(r); }
 153
 154 mp_result mp_rat_denom(mp_rat r, mp_int z) {
 155   return mp_int_copy(MP_DENOM_P(r), z);
 156 }
 157
 158 mp_int mp_rat_denom_ref(mp_rat r) { return MP_DENOM_P(r); }
 159
 160 mp_sign mp_rat_sign(mp_rat r) { return MP_NUMER_SIGN(r); }
 161
 162 mp_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
 173 void 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
 178 mp_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
 189 mp_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
 200 mp_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
 220 mp_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
 224 mp_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
 228 mp_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
 244 mp_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
 282 mp_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
 301 CLEANUP:
 302   mp_int_clear(&tmp);
 303   return res;
 304 }
 305
 306 mp_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
 325 CLEANUP:
 326   mp_int_clear(&tmp);
 327   return res;
 328 }
 329
 330 mp_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
 343 mp_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
 359 mp_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 to
 370      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
 378 int mp_rat_compare(mp_rat a, mp_rat b) {
 379   /* Quick check for opposite signs.  Works because the sign of the numerator
 380      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
 400 int mp_rat_compare_unsigned(mp_rat a, mp_rat b) {
 401   /* If the denominators are equal, we can quickly compare numerators without
 402      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
 432 int mp_rat_compare_zero(mp_rat r) { return mp_int_compare_zero(MP_NUMER_P(r)); }
 433
 434 int 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
 448 CLEANUP:
 449   mp_rat_clear(&tmp);
 450   return out;
 451 }
 452
 453 bool mp_rat_is_integer(mp_rat r) {
 454   return (mp_int_compare_value(MP_DENOM_P(r), 1) == 0);
 455 }
 456
 457 mp_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
 468 mp_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 the
 470      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 exceed
 482      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
 494 mp_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 absolute
 505      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 right
 525      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 is
 542      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 the
 587      displayed digits will be zero due to the precision, a negative shouldn't
 588      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
 621 CLEANUP:
 622   while (--last >= 0) mp_int_clear(TEMP(last));
 623
 624   return res;
 625 }
 626
 627 mp_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 is
 636      included in n_len, and the space for the NUL that is counted in n_len
 637      counts for the separator here.  The space for the NUL counted in d_len
 638      counts for the final terminator here. */
 639
 640   return n_len + d_len;
 641 }
 642
 643 mp_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
 656 mp_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
 660 mp_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 supplied
 697    as zero to use "default" behaviour.
 698
 699    This function will accept either a/b notation or decimal notation.
 700  */
 701 mp_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
 716 mp_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
 720 mp_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 because
 751      mp_int_read_string() will consider whitespace or sign flags to be valid
 752      starting characters for a value, and we do not want them following the
 753      decimal point.
 754
 755      Once we have done this check, it is safe to read in the value of the
 756      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 the
 798        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 the
 820        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 fractional
 837        part was truncated or not.  That information is saved from when we
 838        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 format
 851    and validity of inputs. */
 852
 853 static 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 greater
 863      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
 887 CLEANUP:
 888   mp_int_clear(&gcd);
 889
 890   return res;
 891 }
 892
 893 static 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 */