commit 3c00244d3d0cfa36ad4e6e4f7f6ddb420e905892
parent e7abf41b23140c1ea3c6e00b32f8818b93f9070f
Author: Andres Navarro <canavarro82@gmail.com>
Date: Sat, 23 Apr 2011 09:56:50 -0300
Added rational source code from imath. Still no modification done.
Diffstat:
| A | src/imrat.c | | | 1092 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | src/imrat.h | | | 123 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
2 files changed, 1215 insertions(+), 0 deletions(-)
diff --git a/src/imrat.c b/src/imrat.c
@@ -0,0 +1,1092 @@
+/*
+ Name: imrat.c
+ Purpose: Arbitrary precision rational arithmetic routines.
+ Author: M. J. Fromberger <http://spinning-yarns.org/michael/>
+ Info: $Id: imrat.c 635 2008-01-08 18:19:40Z sting $
+
+ Copyright (C) 2002-2007 Michael J. Fromberger, All Rights Reserved.
+
+ Permission is hereby granted, free of charge, to any person
+ obtaining a copy of this software and associated documentation files
+ (the "Software"), to deal in the Software without restriction,
+ including without limitation the rights to use, copy, modify, merge,
+ publish, distribute, sublicense, and/or sell copies of the Software,
+ and to permit persons to whom the Software is furnished to do so,
+ subject to the following conditions:
+
+ The above copyright notice and this permission notice shall be
+ included in all copies or substantial portions of the Software.
+
+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+ EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+ NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
+ BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
+ ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+ CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ SOFTWARE.
+ */
+
+#include "imrat.h"
+#include <stdlib.h>
+#include <string.h>
+#include <ctype.h>
+#include <assert.h>
+
+/* {{{ Useful macros */
+
+#define TEMP(K) (temp + (K))
+#define SETUP(E, C) \
+do{if((res = (E)) != MP_OK) goto CLEANUP; ++(C);}while(0)
+
+/* Argument checking:
+ Use CHECK() where a return value is required; NRCHECK() elsewhere */
+#define CHECK(TEST) assert(TEST)
+#define NRCHECK(TEST) assert(TEST)
+
+/* }}} */
+
+/* Reduce the given rational, in place, to lowest terms and canonical
+ form. Zero is represented as 0/1, one as 1/1. Signs are adjusted
+ so that the sign of the numerator is definitive. */
+static mp_result s_rat_reduce(mp_rat r);
+
+/* Common code for addition and subtraction operations on rationals. */
+static mp_result s_rat_combine(mp_rat a, mp_rat b, mp_rat c,
+ mp_result (*comb_f)(mp_int, mp_int, mp_int));
+
+/* {{{ mp_rat_init(r) */
+
+mp_result mp_rat_init(mp_rat r)
+{
+ return mp_rat_init_size(r, 0, 0);
+}
+
+/* }}} */
+
+/* {{{ mp_rat_alloc() */
+
+mp_rat mp_rat_alloc(void)
+{
+ mp_rat out = malloc(sizeof(*out));
+
+ if(out != NULL) {
+ if(mp_rat_init(out) != MP_OK) {
+ free(out);
+ return NULL;
+ }
+ }
+
+ return out;
+}
+
+/* }}} */
+
+/* {{{ mp_rat_init_size(r, n_prec, d_prec) */
+
+mp_result mp_rat_init_size(mp_rat r, mp_size n_prec, mp_size d_prec)
+{
+ mp_result res;
+
+ if((res = mp_int_init_size(MP_NUMER_P(r), n_prec)) != MP_OK)
+ return res;
+ if((res = mp_int_init_size(MP_DENOM_P(r), d_prec)) != MP_OK) {
+ mp_int_clear(MP_NUMER_P(r));
+ return res;
+ }
+
+ return mp_int_set_value(MP_DENOM_P(r), 1);
+}
+
+/* }}} */
+
+/* {{{ mp_rat_init_copy(r, old) */
+
+mp_result mp_rat_init_copy(mp_rat r, mp_rat old)
+{
+ mp_result res;
+
+ if((res = mp_int_init_copy(MP_NUMER_P(r), MP_NUMER_P(old))) != MP_OK)
+ return res;
+ if((res = mp_int_init_copy(MP_DENOM_P(r), MP_DENOM_P(old))) != MP_OK)
+ mp_int_clear(MP_NUMER_P(r));
+
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_rat_set_value(r, numer, denom) */
+
+mp_result mp_rat_set_value(mp_rat r, int numer, int denom)
+{
+ mp_result res;
+
+ if(denom == 0)
+ return MP_UNDEF;
+
+ if((res = mp_int_set_value(MP_NUMER_P(r), numer)) != MP_OK)
+ return res;
+ if((res = mp_int_set_value(MP_DENOM_P(r), denom)) != MP_OK)
+ return res;
+
+ return s_rat_reduce(r);
+}
+
+/* }}} */
+
+/* {{{ mp_rat_clear(r) */
+
+void mp_rat_clear(mp_rat r)
+{
+ mp_int_clear(MP_NUMER_P(r));
+ mp_int_clear(MP_DENOM_P(r));
+
+}
+
+/* }}} */
+
+/* {{{ mp_rat_free(r) */
+
+void mp_rat_free(mp_rat r)
+{
+ NRCHECK(r != NULL);
+
+ if(r->num.digits != NULL)
+ mp_rat_clear(r);
+
+ free(r);
+}
+
+/* }}} */
+
+/* {{{ mp_rat_numer(r, z) */
+
+mp_result mp_rat_numer(mp_rat r, mp_int z)
+{
+ return mp_int_copy(MP_NUMER_P(r), z);
+}
+
+/* }}} */
+
+/* {{{ mp_rat_denom(r, z) */
+
+mp_result mp_rat_denom(mp_rat r, mp_int z)
+{
+ return mp_int_copy(MP_DENOM_P(r), z);
+}
+
+/* }}} */
+
+/* {{{ mp_rat_sign(r) */
+
+mp_sign mp_rat_sign(mp_rat r)
+{
+ return MP_SIGN(MP_NUMER_P(r));
+}
+
+/* }}} */
+
+/* {{{ mp_rat_copy(a, c) */
+
+mp_result mp_rat_copy(mp_rat a, mp_rat c)
+{
+ mp_result res;
+
+ if((res = mp_int_copy(MP_NUMER_P(a), MP_NUMER_P(c))) != MP_OK)
+ return res;
+
+ res = mp_int_copy(MP_DENOM_P(a), MP_DENOM_P(c));
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_rat_zero(r) */
+
+void mp_rat_zero(mp_rat r)
+{
+ mp_int_zero(MP_NUMER_P(r));
+ mp_int_set_value(MP_DENOM_P(r), 1);
+
+}
+
+/* }}} */
+
+/* {{{ mp_rat_abs(a, c) */
+
+mp_result mp_rat_abs(mp_rat a, mp_rat c)
+{
+ mp_result res;
+
+ if((res = mp_int_abs(MP_NUMER_P(a), MP_NUMER_P(c))) != MP_OK)
+ return res;
+
+ res = mp_int_abs(MP_DENOM_P(a), MP_DENOM_P(c));
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_rat_neg(a, c) */
+
+mp_result mp_rat_neg(mp_rat a, mp_rat c)
+{
+ mp_result res;
+
+ if((res = mp_int_neg(MP_NUMER_P(a), MP_NUMER_P(c))) != MP_OK)
+ return res;
+
+ res = mp_int_copy(MP_DENOM_P(a), MP_DENOM_P(c));
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_rat_recip(a, c) */
+
+mp_result mp_rat_recip(mp_rat a, mp_rat c)
+{
+ mp_result res;
+
+ if(mp_rat_compare_zero(a) == 0)
+ return MP_UNDEF;
+
+ if((res = mp_rat_copy(a, c)) != MP_OK)
+ return res;
+
+ mp_int_swap(MP_NUMER_P(c), MP_DENOM_P(c));
+
+ /* Restore the signs of the swapped elements */
+ {
+ mp_sign tmp = MP_SIGN(MP_NUMER_P(c));
+
+ MP_SIGN(MP_NUMER_P(c)) = MP_SIGN(MP_DENOM_P(c));
+ MP_SIGN(MP_DENOM_P(c)) = tmp;
+ }
+
+ return MP_OK;
+}
+
+/* }}} */
+
+/* {{{ mp_rat_add(a, b, c) */
+
+mp_result mp_rat_add(mp_rat a, mp_rat b, mp_rat c)
+{
+ return s_rat_combine(a, b, c, mp_int_add);
+
+}
+
+/* }}} */
+
+/* {{{ mp_rat_sub(a, b, c) */
+
+mp_result mp_rat_sub(mp_rat a, mp_rat b, mp_rat c)
+{
+ return s_rat_combine(a, b, c, mp_int_sub);
+
+}
+
+/* }}} */
+
+/* {{{ mp_rat_mul(a, b, c) */
+
+mp_result mp_rat_mul(mp_rat a, mp_rat b, mp_rat c)
+{
+ mp_result res;
+
+ if((res = mp_int_mul(MP_NUMER_P(a), MP_NUMER_P(b), MP_NUMER_P(c))) != MP_OK)
+ return res;
+
+ if(mp_int_compare_zero(MP_NUMER_P(c)) != 0) {
+ if((res = mp_int_mul(MP_DENOM_P(a), MP_DENOM_P(b), MP_DENOM_P(c))) != MP_OK)
+ return res;
+ }
+
+ return s_rat_reduce(c);
+}
+
+/* }}} */
+
+/* {{{ mp_int_div(a, b, c) */
+
+mp_result mp_rat_div(mp_rat a, mp_rat b, mp_rat c)
+{
+ mp_result res = MP_OK;
+
+ if(mp_rat_compare_zero(b) == 0)
+ return MP_UNDEF;
+
+ if(c == a || c == b) {
+ mpz_t tmp;
+
+ if((res = mp_int_init(&tmp)) != MP_OK) return res;
+ if((res = mp_int_mul(MP_NUMER_P(a), MP_DENOM_P(b), &tmp)) != MP_OK)
+ goto CLEANUP;
+ if((res = mp_int_mul(MP_DENOM_P(a), MP_NUMER_P(b), MP_DENOM_P(c))) != MP_OK)
+ goto CLEANUP;
+ res = mp_int_copy(&tmp, MP_NUMER_P(c));
+
+ CLEANUP:
+ mp_int_clear(&tmp);
+ }
+ else {
+ if((res = mp_int_mul(MP_NUMER_P(a), MP_DENOM_P(b), MP_NUMER_P(c))) != MP_OK)
+ return res;
+ if((res = mp_int_mul(MP_DENOM_P(a), MP_NUMER_P(b), MP_DENOM_P(c))) != MP_OK)
+ return res;
+ }
+
+ if(res != MP_OK)
+ return res;
+ else
+ return s_rat_reduce(c);
+}
+
+/* }}} */
+
+/* {{{ mp_rat_add_int(a, b, c) */
+
+mp_result mp_rat_add_int(mp_rat a, mp_int b, mp_rat c)
+{
+ mpz_t tmp;
+ mp_result res;
+
+ if((res = mp_int_init_copy(&tmp, b)) != MP_OK)
+ return res;
+
+ if((res = mp_int_mul(&tmp, MP_DENOM_P(a), &tmp)) != MP_OK)
+ goto CLEANUP;
+
+ if((res = mp_rat_copy(a, c)) != MP_OK)
+ goto CLEANUP;
+
+ if((res = mp_int_add(MP_NUMER_P(c), &tmp, MP_NUMER_P(c))) != MP_OK)
+ goto CLEANUP;
+
+ res = s_rat_reduce(c);
+
+ CLEANUP:
+ mp_int_clear(&tmp);
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_rat_sub_int(a, b, c) */
+
+mp_result mp_rat_sub_int(mp_rat a, mp_int b, mp_rat c)
+{
+ mpz_t tmp;
+ mp_result res;
+
+ if((res = mp_int_init_copy(&tmp, b)) != MP_OK)
+ return res;
+
+ if((res = mp_int_mul(&tmp, MP_DENOM_P(a), &tmp)) != MP_OK)
+ goto CLEANUP;
+
+ if((res = mp_rat_copy(a, c)) != MP_OK)
+ goto CLEANUP;
+
+ if((res = mp_int_sub(MP_NUMER_P(c), &tmp, MP_NUMER_P(c))) != MP_OK)
+ goto CLEANUP;
+
+ res = s_rat_reduce(c);
+
+ CLEANUP:
+ mp_int_clear(&tmp);
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_rat_mul_int(a, b, c) */
+
+mp_result mp_rat_mul_int(mp_rat a, mp_int b, mp_rat c)
+{
+ mp_result res;
+
+ if((res = mp_rat_copy(a, c)) != MP_OK)
+ return res;
+
+ if((res = mp_int_mul(MP_NUMER_P(c), b, MP_NUMER_P(c))) != MP_OK)
+ return res;
+
+ return s_rat_reduce(c);
+}
+
+/* }}} */
+
+/* {{{ mp_rat_div_int(a, b, c) */
+
+mp_result mp_rat_div_int(mp_rat a, mp_int b, mp_rat c)
+{
+ mp_result res;
+
+ if(mp_int_compare_zero(b) == 0)
+ return MP_UNDEF;
+
+ if((res = mp_rat_copy(a, c)) != MP_OK)
+ return res;
+
+ if((res = mp_int_mul(MP_DENOM_P(c), b, MP_DENOM_P(c))) != MP_OK)
+ return res;
+
+ return s_rat_reduce(c);
+}
+
+/* }}} */
+
+/* {{{ mp_rat_expt(a, b, c) */
+
+mp_result mp_rat_expt(mp_rat a, mp_small b, mp_rat c)
+{
+ mp_result res;
+
+ /* Special cases for easy powers. */
+ if(b == 0)
+ return mp_rat_set_value(c, 1, 1);
+ else if(b == 1)
+ return mp_rat_copy(a, c);
+
+ /* Since rationals are always stored in lowest terms, it is not
+ necessary to reduce again when raising to an integer power. */
+ if((res = mp_int_expt(MP_NUMER_P(a), b, MP_NUMER_P(c))) != MP_OK)
+ return res;
+
+ return mp_int_expt(MP_DENOM_P(a), b, MP_DENOM_P(c));
+}
+
+/* }}} */
+
+/* {{{ mp_rat_compare(a, b) */
+
+int mp_rat_compare(mp_rat a, mp_rat b)
+{
+ /* Quick check for opposite signs. Works because the sign of the
+ numerator is always definitive. */
+ if(MP_SIGN(MP_NUMER_P(a)) != MP_SIGN(MP_NUMER_P(b))) {
+ if(MP_SIGN(MP_NUMER_P(a)) == MP_ZPOS)
+ return 1;
+ else
+ return -1;
+ }
+ else {
+ /* Compare absolute magnitudes; if both are positive, the answer
+ stands, otherwise it needs to be reflected about zero. */
+ int cmp = mp_rat_compare_unsigned(a, b);
+
+ if(MP_SIGN(MP_NUMER_P(a)) == MP_ZPOS)
+ return cmp;
+ else
+ return -cmp;
+ }
+}
+
+/* }}} */
+
+/* {{{ mp_rat_compare_unsigned(a, b) */
+
+int mp_rat_compare_unsigned(mp_rat a, mp_rat b)
+{
+ /* If the denominators are equal, we can quickly compare numerators
+ without multiplying. Otherwise, we actually have to do some work. */
+ if(mp_int_compare_unsigned(MP_DENOM_P(a), MP_DENOM_P(b)) == 0)
+ return mp_int_compare_unsigned(MP_NUMER_P(a), MP_NUMER_P(b));
+
+ else {
+ mpz_t temp[2];
+ mp_result res;
+ int cmp = INT_MAX, last = 0;
+
+ /* t0 = num(a) * den(b), t1 = num(b) * den(a) */
+ SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(a)), last);
+ SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(b)), last);
+
+ if((res = mp_int_mul(TEMP(0), MP_DENOM_P(b), TEMP(0))) != MP_OK ||
+ (res = mp_int_mul(TEMP(1), MP_DENOM_P(a), TEMP(1))) != MP_OK)
+ goto CLEANUP;
+
+ cmp = mp_int_compare_unsigned(TEMP(0), TEMP(1));
+
+ CLEANUP:
+ while(--last >= 0)
+ mp_int_clear(TEMP(last));
+
+ return cmp;
+ }
+}
+
+/* }}} */
+
+/* {{{ mp_rat_compare_zero(r) */
+
+int mp_rat_compare_zero(mp_rat r)
+{
+ return mp_int_compare_zero(MP_NUMER_P(r));
+
+}
+
+/* }}} */
+
+/* {{{ mp_rat_compare_value(r, n, d) */
+
+int mp_rat_compare_value(mp_rat r, mp_small n, mp_small d)
+{
+ mpq_t tmp;
+ mp_result res;
+ int out = INT_MAX;
+
+ if((res = mp_rat_init(&tmp)) != MP_OK)
+ return out;
+ if((res = mp_rat_set_value(&tmp, n, d)) != MP_OK)
+ goto CLEANUP;
+
+ out = mp_rat_compare(r, &tmp);
+
+ CLEANUP:
+ mp_rat_clear(&tmp);
+ return out;
+}
+
+/* }}} */
+
+/* {{{ mp_rat_is_integer(r) */
+
+int mp_rat_is_integer(mp_rat r)
+{
+ return (mp_int_compare_value(MP_DENOM_P(r), 1) == 0);
+}
+
+/* }}} */
+
+/* {{{ mp_rat_to_ints(r, *num, *den) */
+
+mp_result mp_rat_to_ints(mp_rat r, mp_small *num, mp_small *den)
+{
+ mp_result res;
+
+ if((res = mp_int_to_int(MP_NUMER_P(r), num)) != MP_OK)
+ return res;
+
+ res = mp_int_to_int(MP_DENOM_P(r), den);
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_rat_to_string(r, radix, *str, limit) */
+
+mp_result mp_rat_to_string(mp_rat r, mp_size radix, char *str, int limit)
+{
+ char *start;
+ int len;
+ mp_result res;
+
+ /* Write the numerator. The sign of the rational number is written
+ by the underlying integer implementation. */
+ if((res = mp_int_to_string(MP_NUMER_P(r), radix, str, limit)) != MP_OK)
+ return res;
+
+ /* If the value is zero, don't bother writing any denominator */
+ if(mp_int_compare_zero(MP_NUMER_P(r)) == 0)
+ return MP_OK;
+
+ /* Locate the end of the numerator, and make sure we are not going to
+ exceed the limit by writing a slash. */
+ len = strlen(str);
+ start = str + len;
+ limit -= len;
+ if(limit == 0)
+ return MP_TRUNC;
+
+ *start++ = '/';
+ limit -= 1;
+
+ res = mp_int_to_string(MP_DENOM_P(r), radix, start, limit);
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_rat_to_decimal(r, radix, prec, *str, limit) */
+mp_result mp_rat_to_decimal(mp_rat r, mp_size radix, mp_size prec,
+ mp_round_mode round, char *str, int limit)
+{
+ mpz_t temp[3];
+ mp_result res;
+ char *start = str;
+ int len, lead_0, left = limit, last = 0;
+
+ SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(r)), last);
+ SETUP(mp_int_init(TEMP(last)), last);
+ SETUP(mp_int_init(TEMP(last)), last);
+
+ /* Get the unsigned integer part by dividing denominator into the
+ absolute value of the numerator. */
+ mp_int_abs(TEMP(0), TEMP(0));
+ if((res = mp_int_div(TEMP(0), MP_DENOM_P(r), TEMP(0), TEMP(1))) != MP_OK)
+ goto CLEANUP;
+
+ /* Now: T0 = integer portion, unsigned;
+ T1 = remainder, from which fractional part is computed. */
+
+ /* Count up leading zeroes after the radix point. */
+ for(lead_0 = 0; lead_0 < prec && mp_int_compare(TEMP(1), MP_DENOM_P(r)) < 0;
+ ++lead_0) {
+ if((res = mp_int_mul_value(TEMP(1), radix, TEMP(1))) != MP_OK)
+ goto CLEANUP;
+ }
+
+ /* Multiply remainder by a power of the radix sufficient to get the
+ right number of significant figures. */
+ if(prec > lead_0) {
+ if((res = mp_int_expt_value(radix, prec - lead_0, TEMP(2))) != MP_OK)
+ goto CLEANUP;
+ if((res = mp_int_mul(TEMP(1), TEMP(2), TEMP(1))) != MP_OK)
+ goto CLEANUP;
+ }
+ if((res = mp_int_div(TEMP(1), MP_DENOM_P(r), TEMP(1), TEMP(2))) != MP_OK)
+ goto CLEANUP;
+
+ /* Now: T1 = significant digits of fractional part;
+ T2 = leftovers, to use for rounding.
+
+ At this point, what we do depends on the rounding mode. The
+ default is MP_ROUND_DOWN, for which everything is as it should be
+ already.
+ */
+ switch(round) {
+ int cmp;
+
+ case MP_ROUND_UP:
+ if(mp_int_compare_zero(TEMP(2)) != 0) {
+ if(prec == 0)
+ res = mp_int_add_value(TEMP(0), 1, TEMP(0));
+ else
+ res = mp_int_add_value(TEMP(1), 1, TEMP(1));
+ }
+ break;
+
+ case MP_ROUND_HALF_UP:
+ case MP_ROUND_HALF_DOWN:
+ if((res = mp_int_mul_pow2(TEMP(2), 1, TEMP(2))) != MP_OK)
+ goto CLEANUP;
+
+ cmp = mp_int_compare(TEMP(2), MP_DENOM_P(r));
+
+ if(round == MP_ROUND_HALF_UP)
+ cmp += 1;
+
+ if(cmp > 0) {
+ if(prec == 0)
+ res = mp_int_add_value(TEMP(0), 1, TEMP(0));
+ else
+ res = mp_int_add_value(TEMP(1), 1, TEMP(1));
+ }
+ break;
+
+ case MP_ROUND_DOWN:
+ break; /* No action required */
+
+ default:
+ return MP_BADARG; /* Invalid rounding specifier */
+ }
+
+ /* The sign of the output should be the sign of the numerator, but
+ if all the displayed digits will be zero due to the precision, a
+ negative shouldn't be shown. */
+ if(MP_SIGN(MP_NUMER_P(r)) == MP_NEG &&
+ (mp_int_compare_zero(TEMP(0)) != 0 ||
+ mp_int_compare_zero(TEMP(1)) != 0)) {
+ *start++ = '-';
+ left -= 1;
+ }
+
+ if((res = mp_int_to_string(TEMP(0), radix, start, left)) != MP_OK)
+ goto CLEANUP;
+
+ len = strlen(start);
+ start += len;
+ left -= len;
+
+ if(prec == 0)
+ goto CLEANUP;
+
+ *start++ = '.';
+ left -= 1;
+
+ if(left < prec + 1) {
+ res = MP_TRUNC;
+ goto CLEANUP;
+ }
+
+ memset(start, '0', lead_0 - 1);
+ left -= lead_0;
+ start += lead_0 - 1;
+
+ res = mp_int_to_string(TEMP(1), radix, start, left);
+
+ CLEANUP:
+ while(--last >= 0)
+ mp_int_clear(TEMP(last));
+
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_rat_string_len(r, radix) */
+
+mp_result mp_rat_string_len(mp_rat r, mp_size radix)
+{
+ mp_result n_len, d_len = 0;
+
+ n_len = mp_int_string_len(MP_NUMER_P(r), radix);
+
+ if(mp_int_compare_zero(MP_NUMER_P(r)) != 0)
+ d_len = mp_int_string_len(MP_DENOM_P(r), radix);
+
+ /* Though simplistic, this formula is correct. Space for the sign
+ flag is included in n_len, and the space for the NUL that is
+ counted in n_len counts for the separator here. The space for
+ the NUL counted in d_len counts for the final terminator here. */
+
+ return n_len + d_len;
+
+}
+
+/* }}} */
+
+/* {{{ mp_rat_decimal_len(r, radix, prec) */
+
+mp_result mp_rat_decimal_len(mp_rat r, mp_size radix, mp_size prec)
+{
+ int z_len, f_len;
+
+ z_len = mp_int_string_len(MP_NUMER_P(r), radix);
+
+ if(prec == 0)
+ f_len = 1; /* terminator only */
+ else
+ f_len = 1 + prec + 1; /* decimal point, digits, terminator */
+
+ return z_len + f_len;
+}
+
+/* }}} */
+
+/* {{{ mp_rat_read_string(r, radix, *str) */
+
+mp_result mp_rat_read_string(mp_rat r, mp_size radix, const char *str)
+{
+ return mp_rat_read_cstring(r, radix, str, NULL);
+}
+
+/* }}} */
+
+/* {{{ mp_rat_read_cstring(r, radix, *str, **end) */
+
+mp_result mp_rat_read_cstring(mp_rat r, mp_size radix, const char *str,
+ char **end)
+{
+ mp_result res;
+ char *endp;
+
+ if((res = mp_int_read_cstring(MP_NUMER_P(r), radix, str, &endp)) != MP_OK &&
+ (res != MP_TRUNC))
+ return res;
+
+ /* Skip whitespace between numerator and (possible) separator */
+ while(isspace((unsigned char) *endp))
+ ++endp;
+
+ /* If there is no separator, we will stop reading at this point. */
+ if(*endp != '/') {
+ mp_int_set_value(MP_DENOM_P(r), 1);
+ if(end != NULL)
+ *end = endp;
+ return res;
+ }
+
+ ++endp; /* skip separator */
+ if((res = mp_int_read_cstring(MP_DENOM_P(r), radix, endp, end)) != MP_OK)
+ return res;
+
+ /* Make sure the value is well-defined */
+ if(mp_int_compare_zero(MP_DENOM_P(r)) == 0)
+ return MP_UNDEF;
+
+ /* Reduce to lowest terms */
+ return s_rat_reduce(r);
+}
+
+/* }}} */
+
+/* {{{ mp_rat_read_ustring(r, radix, *str, **end) */
+
+/* Read a string and figure out what format it's in. The radix may be
+ supplied as zero to use "default" behaviour.
+
+ This function will accept either a/b notation or decimal notation.
+ */
+mp_result mp_rat_read_ustring(mp_rat r, mp_size radix, const char *str,
+ char **end)
+{
+ char *endp;
+ mp_result res;
+
+ if(radix == 0)
+ radix = 10; /* default to decimal input */
+
+ if((res = mp_rat_read_cstring(r, radix, str, &endp)) != MP_OK) {
+ if(res == MP_TRUNC) {
+ if(*endp == '.')
+ res = mp_rat_read_cdecimal(r, radix, str, &endp);
+ }
+ else
+ return res;
+ }
+
+ if(end != NULL)
+ *end = endp;
+
+ return res;
+}
+
+/* }}} */
+
+/* {{{ mp_rat_read_decimal(r, radix, *str) */
+
+mp_result mp_rat_read_decimal(mp_rat r, mp_size radix, const char *str)
+{
+ return mp_rat_read_cdecimal(r, radix, str, NULL);
+}
+
+/* }}} */
+
+/* {{{ mp_rat_read_cdecimal(r, radix, *str, **end) */
+
+mp_result mp_rat_read_cdecimal(mp_rat r, mp_size radix, const char *str,
+ char **end)
+{
+ mp_result res;
+ mp_sign osign;
+ char *endp;
+
+ while(isspace((unsigned char) *str))
+ ++str;
+
+ switch(*str) {
+ case '-':
+ osign = MP_NEG;
+ break;
+ default:
+ osign = MP_ZPOS;
+ }
+
+ if((res = mp_int_read_cstring(MP_NUMER_P(r), radix, str, &endp)) != MP_OK &&
+ (res != MP_TRUNC))
+ return res;
+
+ /* This needs to be here. */
+ (void) mp_int_set_value(MP_DENOM_P(r), 1);
+
+ if(*endp != '.') {
+ if(end != NULL)
+ *end = endp;
+ return res;
+ }
+
+ /* If the character following the decimal point is whitespace or a
+ sign flag, we will consider this a truncated value. This special
+ case is because mp_int_read_string() will consider whitespace or
+ sign flags to be valid starting characters for a value, and we do
+ not want them following the decimal point.
+
+ Once we have done this check, it is safe to read in the value of
+ the fractional piece as a regular old integer.
+ */
+ ++endp;
+ if(*endp == '\0') {
+ if(end != NULL)
+ *end = endp;
+ return MP_OK;
+ }
+ else if(isspace((unsigned char) *endp) || *endp == '-' || *endp == '+') {
+ return MP_TRUNC;
+ }
+ else {
+ mpz_t frac;
+ mp_result save_res;
+ char *save = endp;
+ int num_lz = 0;
+
+ /* Make a temporary to hold the part after the decimal point. */
+ if((res = mp_int_init(&frac)) != MP_OK)
+ return res;
+
+ if((res = mp_int_read_cstring(&frac, radix, endp, &endp)) != MP_OK &&
+ (res != MP_TRUNC))
+ goto CLEANUP;
+
+ /* Save this response for later. */
+ save_res = res;
+
+ if(mp_int_compare_zero(&frac) == 0)
+ goto FINISHED;
+
+ /* Discard trailing zeroes (somewhat inefficiently) */
+ while(mp_int_divisible_value(&frac, radix))
+ if((res = mp_int_div_value(&frac, radix, &frac, NULL)) != MP_OK)
+ goto CLEANUP;
+
+ /* Count leading zeros after the decimal point */
+ while(save[num_lz] == '0')
+ ++num_lz;
+
+ /* Find the least power of the radix that is at least as large as
+ the significant value of the fractional part, ignoring leading
+ zeroes. */
+ (void) mp_int_set_value(MP_DENOM_P(r), radix);
+
+ while(mp_int_compare(MP_DENOM_P(r), &frac) < 0) {
+ if((res = mp_int_mul_value(MP_DENOM_P(r), radix, MP_DENOM_P(r))) != MP_OK)
+ goto CLEANUP;
+ }
+
+ /* Also shift by enough to account for leading zeroes */
+ while(num_lz > 0) {
+ if((res = mp_int_mul_value(MP_DENOM_P(r), radix, MP_DENOM_P(r))) != MP_OK)
+ goto CLEANUP;
+
+ --num_lz;
+ }
+
+ /* Having found this power, shift the numerator leftward that
+ many, digits, and add the nonzero significant digits of the
+ fractional part to get the result. */
+ if((res = mp_int_mul(MP_NUMER_P(r), MP_DENOM_P(r), MP_NUMER_P(r))) != MP_OK)
+ goto CLEANUP;
+
+ { /* This addition needs to be unsigned. */
+ MP_SIGN(MP_NUMER_P(r)) = MP_ZPOS;
+ if((res = mp_int_add(MP_NUMER_P(r), &frac, MP_NUMER_P(r))) != MP_OK)
+ goto CLEANUP;
+
+ MP_SIGN(MP_NUMER_P(r)) = osign;
+ }
+ if((res = s_rat_reduce(r)) != MP_OK)
+ goto CLEANUP;
+
+ /* At this point, what we return depends on whether reading the
+ fractional part was truncated or not. That information is
+ saved from when we called mp_int_read_string() above. */
+ FINISHED:
+ res = save_res;
+ if(end != NULL)
+ *end = endp;
+
+ CLEANUP:
+ mp_int_clear(&frac);
+
+ return res;
+ }
+}
+
+/* }}} */
+
+/* Private functions for internal use. Make unchecked assumptions
+ about format and validity of inputs. */
+
+/* {{{ s_rat_reduce(r) */
+
+static mp_result s_rat_reduce(mp_rat r)
+{
+ mpz_t gcd;
+ mp_result res = MP_OK;
+
+ if(mp_int_compare_zero(MP_NUMER_P(r)) == 0) {
+ mp_int_set_value(MP_DENOM_P(r), 1);
+ return MP_OK;
+ }
+
+ /* If the greatest common divisor of the numerator and denominator
+ is greater than 1, divide it out. */
+ if((res = mp_int_init(&gcd)) != MP_OK)
+ return res;
+
+ if((res = mp_int_gcd(MP_NUMER_P(r), MP_DENOM_P(r), &gcd)) != MP_OK)
+ goto CLEANUP;
+
+ if(mp_int_compare_value(&gcd, 1) != 0) {
+ if((res = mp_int_div(MP_NUMER_P(r), &gcd, MP_NUMER_P(r), NULL)) != MP_OK)
+ goto CLEANUP;
+ if((res = mp_int_div(MP_DENOM_P(r), &gcd, MP_DENOM_P(r), NULL)) != MP_OK)
+ goto CLEANUP;
+ }
+
+ /* Fix up the signs of numerator and denominator */
+ if(MP_SIGN(MP_NUMER_P(r)) == MP_SIGN(MP_DENOM_P(r)))
+ MP_SIGN(MP_NUMER_P(r)) = MP_SIGN(MP_DENOM_P(r)) = MP_ZPOS;
+ else {
+ MP_SIGN(MP_NUMER_P(r)) = MP_NEG;
+ MP_SIGN(MP_DENOM_P(r)) = MP_ZPOS;
+ }
+
+ CLEANUP:
+ mp_int_clear(&gcd);
+
+ return res;
+}
+
+/* }}} */
+
+/* {{{ s_rat_combine(a, b, c, comb_f) */
+
+static mp_result s_rat_combine(mp_rat a, mp_rat b, mp_rat c,
+ mp_result (*comb_f)(mp_int, mp_int, mp_int))
+{
+ mp_result res;
+
+ /* Shortcut when denominators are already common */
+ if(mp_int_compare(MP_DENOM_P(a), MP_DENOM_P(b)) == 0) {
+ if((res = (comb_f)(MP_NUMER_P(a), MP_NUMER_P(b), MP_NUMER_P(c))) != MP_OK)
+ return res;
+ if((res = mp_int_copy(MP_DENOM_P(a), MP_DENOM_P(c))) != MP_OK)
+ return res;
+
+ return s_rat_reduce(c);
+ }
+ else {
+ mpz_t temp[2];
+ int last = 0;
+
+ SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(a)), last);
+ SETUP(mp_int_init_copy(TEMP(last), MP_NUMER_P(b)), last);
+
+ if((res = mp_int_mul(TEMP(0), MP_DENOM_P(b), TEMP(0))) != MP_OK)
+ goto CLEANUP;
+ if((res = mp_int_mul(TEMP(1), MP_DENOM_P(a), TEMP(1))) != MP_OK)
+ goto CLEANUP;
+ if((res = (comb_f)(TEMP(0), TEMP(1), MP_NUMER_P(c))) != MP_OK)
+ goto CLEANUP;
+
+ res = mp_int_mul(MP_DENOM_P(a), MP_DENOM_P(b), MP_DENOM_P(c));
+
+ CLEANUP:
+ while(--last >= 0)
+ mp_int_clear(TEMP(last));
+
+ if(res == MP_OK)
+ return s_rat_reduce(c);
+ else
+ return res;
+ }
+}
+
+/* }}} */
+
+/* Here there be dragons */
diff --git a/src/imrat.h b/src/imrat.h
@@ -0,0 +1,123 @@
+/*
+ Name: imrat.h
+ Purpose: Arbitrary precision rational arithmetic routines.
+ Author: M. J. Fromberger <http://spinning-yarns.org/michael/>
+ Info: $Id: imrat.h 635 2008-01-08 18:19:40Z sting $
+
+ Copyright (C) 2002-2007 Michael J. Fromberger, All Rights Reserved.
+
+ Permission is hereby granted, free of charge, to any person
+ obtaining a copy of this software and associated documentation files
+ (the "Software"), to deal in the Software without restriction,
+ including without limitation the rights to use, copy, modify, merge,
+ publish, distribute, sublicense, and/or sell copies of the Software,
+ and to permit persons to whom the Software is furnished to do so,
+ subject to the following conditions:
+
+ The above copyright notice and this permission notice shall be
+ included in all copies or substantial portions of the Software.
+
+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+ EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+ NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
+ BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
+ ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+ CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ SOFTWARE.
+ */
+
+#ifndef IMRAT_H_
+#define IMRAT_H_
+
+#include "imath.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef struct mpq {
+ mpz_t num; /* Numerator */
+ mpz_t den; /* Denominator, <> 0 */
+} mpq_t, *mp_rat;
+
+#define MP_NUMER_P(Q) (&((Q)->num)) /* Pointer to numerator */
+#define MP_DENOM_P(Q) (&((Q)->den)) /* Pointer to denominator */
+
+/* Rounding constants */
+typedef enum {
+ MP_ROUND_DOWN,
+ MP_ROUND_HALF_UP,
+ MP_ROUND_UP,
+ MP_ROUND_HALF_DOWN
+} mp_round_mode;
+
+mp_result mp_rat_init(mp_rat r);
+mp_rat mp_rat_alloc(void);
+mp_result mp_rat_init_size(mp_rat r, mp_size n_prec, mp_size d_prec);
+mp_result mp_rat_init_copy(mp_rat r, mp_rat old);
+mp_result mp_rat_set_value(mp_rat r, int numer, int denom);
+void mp_rat_clear(mp_rat r);
+void mp_rat_free(mp_rat r);
+mp_result mp_rat_numer(mp_rat r, mp_int z); /* z = num(r) */
+mp_result mp_rat_denom(mp_rat r, mp_int z); /* z = den(r) */
+mp_sign mp_rat_sign(mp_rat r);
+
+mp_result mp_rat_copy(mp_rat a, mp_rat c); /* c = a */
+void mp_rat_zero(mp_rat r); /* r = 0 */
+mp_result mp_rat_abs(mp_rat a, mp_rat c); /* c = |a| */
+mp_result mp_rat_neg(mp_rat a, mp_rat c); /* c = -a */
+mp_result mp_rat_recip(mp_rat a, mp_rat c); /* c = 1 / a */
+mp_result mp_rat_add(mp_rat a, mp_rat b, mp_rat c); /* c = a + b */
+mp_result mp_rat_sub(mp_rat a, mp_rat b, mp_rat c); /* c = a - b */
+mp_result mp_rat_mul(mp_rat a, mp_rat b, mp_rat c); /* c = a * b */
+mp_result mp_rat_div(mp_rat a, mp_rat b, mp_rat c); /* c = a / b */
+
+mp_result mp_rat_add_int(mp_rat a, mp_int b, mp_rat c); /* c = a + b */
+mp_result mp_rat_sub_int(mp_rat a, mp_int b, mp_rat c); /* c = a - b */
+mp_result mp_rat_mul_int(mp_rat a, mp_int b, mp_rat c); /* c = a * b */
+mp_result mp_rat_div_int(mp_rat a, mp_int b, mp_rat c); /* c = a / b */
+mp_result mp_rat_expt(mp_rat a, mp_small b, mp_rat c); /* c = a ^ b */
+
+int mp_rat_compare(mp_rat a, mp_rat b); /* a <=> b */
+int mp_rat_compare_unsigned(mp_rat a, mp_rat b); /* |a| <=> |b| */
+int mp_rat_compare_zero(mp_rat r); /* r <=> 0 */
+int mp_rat_compare_value(mp_rat r, mp_small n, mp_small d); /* r <=> n/d */
+int mp_rat_is_integer(mp_rat r);
+
+/* Convert to integers, if representable (returns MP_RANGE if not). */
+mp_result mp_rat_to_ints(mp_rat r, mp_small *num, mp_small *den);
+
+/* Convert to nul-terminated string with the specified radix, writing
+ at most limit characters including the nul terminator. */
+mp_result mp_rat_to_string(mp_rat r, mp_size radix, char *str, int limit);
+
+/* Convert to decimal format in the specified radix and precision,
+ writing at most limit characters including a nul terminator. */
+mp_result mp_rat_to_decimal(mp_rat r, mp_size radix, mp_size prec,
+ mp_round_mode round, char *str, int limit);
+
+/* Return the number of characters required to represent r in the given
+ radix. May over-estimate. */
+mp_result mp_rat_string_len(mp_rat r, mp_size radix);
+
+/* Return the number of characters required to represent r in decimal
+ format with the given radix and precision. May over-estimate. */
+mp_result mp_rat_decimal_len(mp_rat r, mp_size radix, mp_size prec);
+
+/* Read zero-terminated string into r */
+mp_result mp_rat_read_string(mp_rat r, mp_size radix, const char *str);
+mp_result mp_rat_read_cstring(mp_rat r, mp_size radix, const char *str,
+ char **end);
+mp_result mp_rat_read_ustring(mp_rat r, mp_size radix, const char *str,
+ char **end);
+
+/* Read zero-terminated string in decimal format into r */
+mp_result mp_rat_read_decimal(mp_rat r, mp_size radix, const char *str);
+mp_result mp_rat_read_cdecimal(mp_rat r, mp_size radix, const char *str,
+ char **end);
+
+#ifdef __cplusplus
+}
+#endif
+#endif /* IMRAT_H_ */
