commit 9ef33b000b7370bc31326fb5b889b5ab740ac042
parent 0945c7a9e3849e431c0ae9501b2f10fcaa7d1d83
Author: Andres Navarro <canavarro82@gmail.com>
Date: Tue, 26 Apr 2011 11:43:11 -0300
Added rationalize & simplest rational to the ground environment.
Diffstat:
5 files changed, 331 insertions(+), 9 deletions(-)
diff --git a/src/kgnumbers.c b/src/kgnumbers.c
@@ -442,7 +442,75 @@ TValue knum_real_to_integer(klisp_State *K, TValue n, kround_mode mode)
klispE_throw(K, "round: infinite value");
return KINERT;
default:
- klispE_throw(K, "denominator: unsopported type");
+ klispE_throw(K, "round: unsopported type");
+ return KINERT;
+ }
+}
+
+TValue knum_simplest_rational(klisp_State *K, TValue n1, TValue n2)
+{
+ /* first check that case that n1 > n2 */
+ if (knum_gtp(K, n1, n2)) {
+ klispE_throw(K, "simplest_rational: result with no primary value "
+ "(n1 > n2)");
+ return KINERT;
+ }
+
+ /* we know that n1 <= n2 */
+ switch(max_ttype(n1, n2)) {
+ case K_TFIXINT:
+ case K_TBIGINT: /* for now do all with bigrat */
+ case K_TBIGRAT: {
+ /* we know that n1 <= n2 */
+ kensure_bigrat(n1);
+ kensure_bigrat(n2);
+ return kbigrat_simplest_rational(K, n1, n2);
+ }
+ case K_TEINF:
+ /* we know that n1 <= n2 */
+ if (tv_equal(n1, n2)) {
+ klispE_throw(K, "simplest rational: result with no primary value");
+ return KINERT;
+ } else if (knegativep(n1) && kpositivep(n2)) {
+ return i2tv(0);
+ } else if (knegativep(n1)) {
+ /* n1 -inf, n2 finite negative */
+ /* ASK John: is this behaviour for infinities ok? */
+ /* Also in the report example both 1/3 & 1/2 are simpler than
+ 2/5... */
+ return knum_real_to_integer(K, n2, K_FLOOR);
+ } else {
+ /* n1 finite positive, n2 +inf */
+ /* ASK John: is this behaviour for infinities ok? */
+ return knum_real_to_integer(K, n1, K_CEILING);
+ }
+ default:
+ klispE_throw(K, "simplest rational: unsopported type");
+ return KINERT;
+ }
+}
+
+TValue knum_rationalize(klisp_State *K, TValue n1, TValue n2)
+{
+ switch(max_ttype(n1, n2)) {
+ case K_TFIXINT:
+ case K_TBIGINT: /* for now do all with bigrat */
+ case K_TBIGRAT: {
+ /* we know that n1 <= n2 */
+ kensure_bigrat(n1);
+ kensure_bigrat(n2);
+ return kbigrat_rationalize(K, n1, n2);
+ }
+ case K_TEINF:
+ if (kfinitep(n1) || !kfinitep(n2)) {
+ return i2tv(0);
+ } else { /* infinite n1, finite n2 */
+ /* ASK John: is this behaviour for infinities ok? */
+ klispE_throw(K, "rationalize: result with no primary value");
+ return KINERT;
+ }
+ default:
+ klispE_throw(K, "rationalize: unsopported type");
return KINERT;
}
}
@@ -1142,4 +1210,28 @@ void kreal_to_integer(klisp_State *K, TValue *xparams, TValue ptree,
}
/* 12.8.5 rationalize, simplest-rational */
-/* TODO */
+void krationalize(klisp_State *K, TValue *xparams, TValue ptree,
+ TValue denv)
+{
+ UNUSED(denv);
+ UNUSED(xparams);
+
+ bind_2tp(K, "rationalize", ptree, "real", krealp, n1,
+ "real", krealp, n2);
+
+ TValue res = knum_rationalize(K, n1, n2);
+ kapply_cc(K, res);
+}
+
+void ksimplest_rational(klisp_State *K, TValue *xparams, TValue ptree,
+ TValue denv)
+{
+ UNUSED(denv);
+ UNUSED(xparams);
+
+ bind_2tp(K, "simplest-rational", ptree, "real", krealp, n1,
+ "real", krealp, n2);
+
+ TValue res = knum_simplest_rational(K, n1, n2);
+ kapply_cc(K, res);
+}
diff --git a/src/kgnumbers.h b/src/kgnumbers.h
@@ -126,7 +126,11 @@ void kreal_to_integer(klisp_State *K, TValue *xparams, TValue ptree,
TValue denv);
/* 12.8.5 rationalize, simplest-rational */
-/* TODO */
+void krationalize(klisp_State *K, TValue *xparams, TValue ptree,
+ TValue denv);
+
+void ksimplest_rational(klisp_State *K, TValue *xparams, TValue ptree,
+ TValue denv);
/* REFACTOR: These should be in a knumber.h header */
diff --git a/src/kground.c b/src/kground.c
@@ -757,7 +757,8 @@ void kinit_ground_env(klisp_State *K)
symbol, i2tv((int32_t) K_ROUND_EVEN));
/* 12.8.5 rationalize, simplest-rational */
- /* TODO */
+ add_applicative(K, ground_env, "rationalize", krationalize, 0);
+ add_applicative(K, ground_env, "simplest-rational", ksimplest_rational, 0);
/*
**
diff --git a/src/krational.c b/src/krational.c
@@ -300,3 +300,231 @@ TValue kbigrat_to_integer(klisp_State *K, TValue tv_bigrat, kround_mode mode)
krooted_tvs_pop(K);
return kbigint_try_fixint(K, tv_quot);
}
+
+/*
+** SOURCE NOTE: this implementation is from the Haskell 98 report
+*/
+/*
+approxRational x eps = simplest (x-eps) (x+eps)
+ where simplest x y | y < x = simplest y x
+ | x == y = xr
+ | x > 0 = simplest' n d n' d'
+ | y < 0 = - simplest' (-n') d' (-n) d
+ | otherwise = 0 :% 1
+ where xr@(n:%d) = toRational x
+ (n':%d') = toRational y
+
+ simplest' n d n' d' -- assumes 0 < n%d < n'%d'
+ | r == 0 = q :% 1
+ | q /= q' = (q+1) :% 1
+ | otherwise = (q*n''+d'') :% n''
+ where (q,r) = quotRem n d
+ (q',r') = quotRem n' d'
+ (n'':%d'') = simplest' d' r' d r
+
+*/
+
+/*
+** NOTE: this reads almost like a Haskell commercial.
+** The c code is an order of magnitude longer. Some of this has to do
+** with the representation of values, some because this is iterative,
+** some because of GC rooting, some because of lack of powerful pattern
+** matching, and so on, and so on
+*/
+
+/* Assumes 0 < n1/d1 < n2/d2 */
+/* GC: Assumes n1, d1, n2, d2, and res are fresh (can be mutated) and rooted */
+static void simplest(klisp_State *K, TValue tv_n1, TValue tv_d1,
+ TValue tv_n2, TValue tv_d2, TValue tv_res)
+{
+ Bigint *n1 = tv2bigint(tv_n1);
+ Bigint *d1 = tv2bigint(tv_d1);
+ Bigint *n2 = tv2bigint(tv_n2);
+ Bigint *d2 = tv2bigint(tv_d2);
+
+ Bigrat *res = tv2bigrat(tv_res);
+
+ /* create vars q1, r1, q2 & r2 */
+ TValue tv_q1 = kbigint_make_simple(K);
+ krooted_tvs_push(K, tv_q1);
+ Bigint *q1 = tv2bigint(tv_q1);
+
+ TValue tv_r1 = kbigint_make_simple(K);
+ krooted_tvs_push(K, tv_r1);
+ Bigint *r1 = tv2bigint(tv_r1);
+
+ TValue tv_q2 = kbigint_make_simple(K);
+ krooted_tvs_push(K, tv_q2);
+ Bigint *q2 = tv2bigint(tv_q2);
+
+ TValue tv_r2 = kbigint_make_simple(K);
+ krooted_tvs_push(K, tv_r2);
+ Bigint *r2 = tv2bigint(tv_r2);
+
+ while(true) {
+ UNUSED(mp_int_div(K, n1, d1, q1, r1));
+ UNUSED(mp_int_div(K, n2, d2, q2, r2));
+
+ if (mp_int_compare_zero(r1) == 0) {
+ /* res = q1 / 1 */
+ mp_rat_zero(K, res);
+ mp_rat_add_int(K, res, q1, res);
+ break;
+ } else if (mp_int_compare(q1, q2) != 0) {
+ /* res = (q1 + 1) / 1 */
+ mp_rat_zero(K, res);
+ mp_int_add_value(K, q1, 1, q1);
+ mp_rat_add_int(K, res, q1, res);
+ break;
+ } else {
+ /* simulate a recursive call */
+ TValue saved_q1 = kbigint_make_simple(K);
+ krooted_tvs_push(K, saved_q1);
+ UNUSED(mp_int_copy(K, q1, tv2bigint(saved_q1)));
+ ks_spush(K, saved_q1);
+ krooted_tvs_pop(K);
+
+ UNUSED(mp_int_copy(K, d2, n1));
+ UNUSED(mp_int_copy(K, d1, n2));
+ UNUSED(mp_int_copy(K, r2, d1));
+ UNUSED(mp_int_copy(K, r1, d2));
+ } /* fall through */
+ }
+
+ /* now, if there were "recursive" calls, complete them */
+ while(!ks_sisempty(K)) {
+ TValue saved_q1 = ks_sget(K);
+ TValue tv_tmp = kbigrat_make_simple(K);
+ krooted_tvs_push(K, tv_tmp);
+ Bigrat *tmp = tv2bigrat(tv_tmp);
+
+ UNUSED(mp_rat_copy(K, res, tmp));
+ /* res = (saved_q * n(res)) + d(res)) / n(res) */
+ UNUSED(mp_rat_zero(K, res));
+ UNUSED(mp_rat_add_int(K, res, tv2bigint(saved_q1), res));
+ UNUSED(mp_rat_mul_int(K, res, MP_NUMER_P(tmp), res));
+ UNUSED(mp_rat_add_int(K, res, MP_DENOM_P(tmp), res));
+ UNUSED(mp_rat_div_int(K, res, MP_NUMER_P(tmp), res));
+ krooted_tvs_pop(K); /* tmp */
+ ks_sdpop(K); /* saved_q1 */
+ }
+
+ krooted_tvs_pop(K); /* q1, r1, q2, r2 */
+ krooted_tvs_pop(K);
+ krooted_tvs_pop(K);
+ krooted_tvs_pop(K);
+
+ return;
+}
+
+TValue kbigrat_simplest_rational(klisp_State *K, TValue tv_n1, TValue tv_n2)
+{
+ TValue tv_res = kbigrat_make_simple(K);
+ krooted_tvs_push(K, tv_res);
+ Bigrat *res = tv2bigrat(tv_res);
+ Bigrat *n1 = tv2bigrat(tv_n1);
+ Bigrat *n2 = tv2bigrat(tv_n2);
+
+ int32_t cmp = mp_rat_compare(K, n1, n2);
+ if (cmp > 0) { /* n1 > n2, swap */
+ TValue temp = tv_n1;
+ tv_n1 = tv_n2;
+ tv_n2 = temp;
+ n1 = tv2bigrat(tv_n1);
+ n2 = tv2bigrat(tv_n2);
+ /* fall through */
+ } else if (cmp == 0) { /* n1 == n2 */
+ krooted_tvs_pop(K);
+ return tv_n1;
+ } /* else fall through */
+
+ /* we now know that n1 < n2 */
+ if (mp_rat_compare_zero(n1) > 0) {
+ /* 0 > n1 > n2 */
+ TValue num1 = kbigint_make_simple(K);
+ krooted_tvs_push(K, num1);
+ UNUSED(mp_rat_numer(K, n1, tv2bigint(num1)));
+
+ TValue den1 = kbigint_make_simple(K);
+ krooted_tvs_push(K, den1);
+ UNUSED(mp_rat_denom(K, n1, tv2bigint(den1)));
+
+ TValue num2 = kbigint_make_simple(K);
+ krooted_tvs_push(K, num2);
+ UNUSED(mp_rat_numer(K, n2, tv2bigint(num2)));
+
+ TValue den2 = kbigint_make_simple(K);
+ krooted_tvs_push(K, den2);
+ UNUSED(mp_rat_denom(K, n2, tv2bigint(den2)));
+
+ simplest(K, num1, den1, num2, den2, tv_res);
+
+ krooted_tvs_pop(K); /* num1, num2, den1, den2 */
+ krooted_tvs_pop(K);
+ krooted_tvs_pop(K);
+ krooted_tvs_pop(K);
+
+ krooted_tvs_pop(K); /* tv_res */
+ return kbigrat_try_integer(K, tv_res);
+ } else if (mp_rat_compare_zero(n2) < 0) {
+ /* n1 < n2 < 0 */
+
+ /* do -(simplest -n2/d2 -n1/d1) */
+
+ TValue num1 = kbigint_make_simple(K);
+ krooted_tvs_push(K, num1);
+ UNUSED(mp_int_neg(K, MP_NUMER_P(n2), tv2bigint(num1)));
+
+ TValue den1 = kbigint_make_simple(K);
+ krooted_tvs_push(K, den1);
+ UNUSED(mp_rat_denom(K, n2, tv2bigint(den1)));
+
+ TValue num2 = kbigint_make_simple(K);
+ krooted_tvs_push(K, num2);
+ UNUSED(mp_int_neg(K, MP_NUMER_P(n1), tv2bigint(num2)));
+
+ TValue den2 = kbigint_make_simple(K);
+ krooted_tvs_push(K, den2);
+ UNUSED(mp_rat_denom(K, n1, tv2bigint(den2)));
+
+ simplest(K, num1, den1, num2, den2, tv_res);
+ mp_rat_neg(K, res, res);
+
+ krooted_tvs_pop(K); /* num1, num2, den1, den2 */
+ krooted_tvs_pop(K);
+ krooted_tvs_pop(K);
+ krooted_tvs_pop(K);
+
+ krooted_tvs_pop(K); /* tv_res */
+ return kbigrat_try_integer(K, tv_res);
+ } else {
+ /* n1 < 0 < n2 */
+ krooted_tvs_pop(K);
+ return i2tv(0);
+ }
+}
+
+TValue kbigrat_rationalize(klisp_State *K, TValue tv_n1, TValue tv_n2)
+{
+ /* delegate all work to simplest_rational */
+ TValue tv_min = kbigrat_make_simple(K);
+ krooted_tvs_push(K, tv_min);
+ TValue tv_max = kbigrat_make_simple(K);
+ krooted_tvs_push(K, tv_max);
+
+ Bigrat *n1 = tv2bigrat(tv_n1);
+ Bigrat *n2 = tv2bigrat(tv_n2);
+ /* it doesn't matter if these are reversed */
+ Bigrat *min = tv2bigrat(tv_min);
+ Bigrat *max = tv2bigrat(tv_max);
+
+ UNUSED(mp_rat_sub(K, n1, n2, min));
+ UNUSED(mp_rat_add(K, n1, n2, max));
+
+ TValue res = kbigrat_simplest_rational(K, tv_min, tv_max);
+
+ krooted_tvs_pop(K);
+ krooted_tvs_pop(K);
+
+ return res;
+}
diff --git a/src/krational.h b/src/krational.h
@@ -166,13 +166,10 @@ TValue kbigrat_abs(klisp_State *K, TValue tv_bigrat);
TValue kbigrat_numerator(klisp_State *K, TValue tv_bigrat);
TValue kbigrat_denominator(klisp_State *K, TValue tv_bigrat);
-/* TODO implement these */
typedef enum { K_FLOOR, K_CEILING, K_TRUNCATE, K_ROUND_EVEN } kround_mode;
TValue kbigrat_to_integer(klisp_State *K, TValue tv_bigrat, kround_mode mode);
-#if 0
-TValue kbigrat_simplest_rational(klisp_State *K, TValue n1, TValue n2);
-TValue kbigrat_rationalize(klisp_State *K, TValue n1, TValue n2);
-#endif
+TValue kbigrat_simplest_rational(klisp_State *K, TValue tv_n1, TValue tv_n2);
+TValue kbigrat_rationalize(klisp_State *K, TValue tv_n1, TValue tv_n2);
#endif
