commit 528609859fad4cefee74f35935cf3a40d8a91a20
parent a795b4251cf9be5a65c1ec000ded006a6c73cbb5
Author: Andres Navarro <canavarro82@gmail.com>
Date: Sat, 23 Apr 2011 10:49:55 -0300
Merged the work I started in the default branch so it is all in this rational branch
Diffstat:
| M | src/imath.c | | | 13 | +++++++------ |
| M | src/imath.h | | | 7 | ++++--- |
| A | src/imrat.c | | | 1081 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | src/imrat.h | | | 139 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| M | src/kobject.h | | | 18 | +++++++++++++++++- |
5 files changed, 1248 insertions(+), 10 deletions(-)
diff --git a/src/imath.c b/src/imath.c
@@ -23,18 +23,19 @@
#include <assert.h>
-#if DEBUG
-#define STATIC /* public */
-#else
-#define STATIC static
-#endif
-
/* Andres Navarro: klisp includes */
#include "kobject.h"
#include "kstate.h"
#include "kmem.h"
#include "kerror.h"
+
+#if DEBUG
+#define STATIC /* public */
+#else
+#define STATIC static
+#endif
+
/* {{{ Constants */
const mp_result MP_OK = 0; /* no error, all is well */
diff --git a/src/imath.h b/src/imath.h
@@ -57,7 +57,9 @@ typedef unsigned int mp_word;
#endif /* USE_C99 */
/* Andres Navarro: Use kobject type instead */
-/*
+typedef Bigint mpz_t, *mp_int;
+
+#if 0
typedef struct mpz {
mp_digit single;
mp_digit *digits;
@@ -65,8 +67,7 @@ typedef struct mpz {
mp_size used;
mp_sign sign;
} mpz_t, *mp_int;
-*/
-typedef Bigint mpz_t, *mp_int;
+#endif
#define MP_SINGLE(Z) ((Z)->single) /* added to correct check in mp_int_clear */
#define MP_DIGITS(Z) ((Z)->digits)
diff --git a/src/imrat.c b/src/imrat.c
@@ -0,0 +1,1081 @@
+/*
+** imrat.c
+** Arbitrary precision rational arithmetic routines.
+** See Copyright Notice in klisp.h
+*/
+
+/*
+** SOURCE NOTE: This is mostly from the IMath library, written by
+** M.J. Fromberger. It is adapted to klisp, mainly in the use of the
+** klisp allocator and fixing of digit size to 32 bits.
+** Imported from version (1.15) updated 01-Feb-2011 at 03:10 PM.
+*/
+
+#include "imrat.h"
+#include <stdlib.h>
+#include <string.h>
+#include <ctype.h>
+#include <assert.h>
+
+/* Andres Navarro: klisp includes */
+#include "kobject.h"
+#include "kstate.h"
+#include "kmem.h"
+#include "kerror.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,139 @@
+/*
+** imrat.h
+** Arbitrary precision rational arithmetic routines.
+** See Copyright Notice in klisp.h
+*/
+
+/*
+** SOURCE NOTE: This is mostly from the IMath library, written by
+** M.J. Fromberger. It is adapted to klisp, mainly in the use of the
+** klisp allocator and fixing of digit size to 32 bits.
+** Imported from version (1.15) updated 01-Feb-2011 at 03:10 PM.
+*/
+
+#ifndef IMRAT_H_
+#define IMRAT_H_
+
+#include "imath.h"
+
+/* Andres Navarro: klisp includes */
+#include "kobject.h"
+#include "kstate.h"
+
+#ifdef USE_C99
+#include <stdint.h>
+#endif
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/* Andres Navarro: Use kobject type instead */
+typedef Bigrat mpq_t, *mp_rat;
+
+#if 0
+typedef struct mpq {
+ mpz_t num; /* Numerator */
+ mpz_t den; /* Denominator, <> 0 */
+} mpq_t, *mp_rat;
+#endif
+
+#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(klisp_State *K, mp_rat r);
+mp_rat mp_rat_alloc(klisp_State *K);
+mp_result mp_rat_init_size(klisp_State *K, mp_rat r, mp_size n_prec,
+ mp_size d_prec);
+mp_result mp_rat_init_copy(klisp_State *K, mp_rat r, mp_rat old);
+mp_result mp_rat_set_value(klisp_State *K, mp_rat r, int numer, int denom);
+void mp_rat_clear(klisp_State *K, mp_rat r);
+void mp_rat_free(klisp_State *K, mp_rat r);
+mp_result mp_rat_numer(klisp_State *K, mp_rat r, mp_int z); /* z = num(r) */
+mp_result mp_rat_denom(klisp_State *K, mp_rat r, mp_int z); /* z = den(r) */
+/* NOTE: this doesn't use the allocator */
+mp_sign mp_rat_sign(mp_rat r);
+
+mp_result mp_rat_copy(klisp_State *K, mp_rat a, mp_rat c); /* c = a */
+/* NOTE: this doesn't use the allocator */
+void mp_rat_zero(mp_rat r); /* r = 0 */
+mp_result mp_rat_abs(klisp_State *K, mp_rat a, mp_rat c); /* c = |a| */
+mp_result mp_rat_neg(klisp_State *K, mp_rat a, mp_rat c); /* c = -a */
+mp_result mp_rat_recip(klisp_State *K, mp_rat a, mp_rat c); /* c = 1 / a */
+/* c = a + b */
+mp_result mp_rat_add(klisp_State *K, mp_rat a, mp_rat b, mp_rat c);
+/* c = a - b */
+mp_result mp_rat_sub(klisp_State *K, mp_rat a, mp_rat b, mp_rat c);
+/* c = a * b */
+mp_result mp_rat_mul(klisp_State *K, mp_rat a, mp_rat b, mp_rat c);
+/* c = a / b */
+mp_result mp_rat_div(klisp_State *K, mp_rat a, mp_rat b, mp_rat c);
+
+/* c = a + b */
+mp_result mp_rat_add_int(klisp_State *K, mp_rat a, mp_int b, mp_rat c);
+/* c = a - b */
+mp_result mp_rat_sub_int(klisp_State *K, mp_rat a, mp_int b, mp_rat c);
+/* c = a * b */
+mp_result mp_rat_mul_int(klisp_State *K, mp_rat a, mp_int b, mp_rat c);
+/* c = a / b */
+mp_result mp_rat_div_int(klisp_State *K, mp_rat a, mp_int b, mp_rat c);
+/* c = a ^ b */
+mp_result mp_rat_expt(klisp_State *K, mp_rat a, mp_small b, mp_rat c);
+
+/* NOTE: because we may need to do multiplications, some of
+ these take a klisp_State */
+int mp_rat_compare(klisp_State *K, mp_rat a, mp_rat b); /* a <=> b */
+/* |a| <=> |b| */
+int mp_rat_compare_unsigned(klisp_State *K, mp_rat a, mp_rat b);
+int mp_rat_compare_zero(mp_rat r); /* r <=> 0 */
+int mp_rat_compare_value(klisp_State *K, 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). */
+/* NOTE: this doesn't use the allocator */
+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(klisp_State *K, mp_rat r, mp_size radix,
+ const char *str);
+mp_result mp_rat_read_cstring(klisp_State *K, mp_rat r, mp_size radix,
+ const char *str, char **end);
+mp_result mp_rat_read_ustring(klisp_State *K, 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(klisp_State *K, mp_rat r, mp_size radix,
+ const char *str);
+mp_result mp_rat_read_cdecimal(klisp_State *K, mp_rat r, mp_size radix,
+ const char *str, char **end);
+
+#ifdef __cplusplus
+}
+#endif
+#endif /* IMRAT_H_ */
diff --git a/src/kobject.h b/src/kobject.h
@@ -179,10 +179,11 @@ typedef struct __attribute__ ((__packed__)) GCheader {
**
** - decide if inexact infinities and reals with no
** primary values are included in K_TDOUBLE
-** - For now we will only use fixints, bigints and exact infinities
+** - For now we will only use fixints, bigints, bigrats and exact infinities
*/
#define K_TAG_FIXINT K_MAKE_VTAG(K_TFIXINT)
#define K_TAG_BIGINT K_MAKE_VTAG(K_TBIGINT)
+#define K_TAG_BIGRAT K_MAKE_VTAG(K_TBIGRAT)
#define K_TAG_EINF K_MAKE_VTAG(K_TEINF)
#define K_TAG_IINF K_MAKE_VTAG(K_TIINF)
@@ -232,6 +233,10 @@ typedef struct __attribute__ ((__packed__)) GCheader {
#define ttisbigint(o) (tbasetype_(o) == K_TAG_BIGINT)
#define ttisinteger(o_) ({ int32_t t_ = tbasetype_(o_); \
t_ == K_TAG_FIXINT || t_ == K_TAG_BIGINT;})
+#define ttisbigrat(o) (tbasetype_(o) == K_TAG_BIGRAT)
+#define ttisrational(o) ({ int32_t t_ = tbasetype_(o_); \
+ t_ == K_TAG_BIGRAT || t_== K_TAG_BIGINT || \
+ t == K_TAG_FIXINT;})
#define ttisnumber(o) (ttype(o) <= K_LAST_NUMBER_TYPE); })
#define ttiseinf(o) (tbasetype_(o) == K_TAG_EINF)
#define ttisiinf(o) (tbasetype_(o) == K_TAG_IINF)
@@ -314,6 +319,17 @@ typedef struct __attribute__ ((__packed__)) {
unsigned char sign;
} Bigint;
+/* NOTE: Notice that both num and den aren't pointers, so, in general, to get
+ the denominator or numerator we have to make a copy, this comes from IMath.
+ If written for klisp I would have put pointers. Maybe I'll later change it
+ but for now minimal ammount of modification to IMath is desired */
+typedef struct __attribute__ ((__packed__)) {
+ CommonHeader;
+/* This is from IMath */
+ Bigint num; /* Numerator */
+ Bigint den; /* Denominator, <> 0 */
+} Bigrat;
+
/* REFACTOR: move these macros somewhere else */
/* NOTE: The use of the intermediate KCONCAT is needed to allow
expansion of the __LINE__ macro. */
