commit 72f839f5a13fd5f05a86517fa4e4ebf205831968
parent cd73337489f1dab7beaf04d13b9b53d1fce79c77
Author: Andres Navarro <canavarro82@gmail.com>
Date: Sun, 15 May 2011 11:27:18 -0300
Added support for inexact reals to div-and-mod & div0-and-mod0.
Diffstat:
3 files changed, 112 insertions(+), 1 deletion(-)
diff --git a/src/kgnumbers.c b/src/kgnumbers.c
@@ -1054,6 +1054,8 @@ void kdiv_mod(klisp_State *K, TValue *xparams, TValue ptree, TValue denv)
TValue tv_div, tv_mod;
+ kensure_same_exactness(K, tv_n, tv_d);
+
if (kfast_zerop(tv_d)) {
klispE_throw_simple(K, "division by zero");
return;
@@ -1096,6 +1098,17 @@ void kdiv_mod(klisp_State *K, TValue *xparams, TValue ptree, TValue denv)
else
tv_div = kbigrat_div0_mod0(K, tv_n, tv_d, &tv_mod);
break;
+ case K_TDOUBLE: {
+ /* both are double */
+ double div, mod;
+ if ((flags & FDIV_ZERO) == 0)
+ div = kdouble_div_mod(dvalue(tv_n), dvalue(tv_d), &mod);
+ else
+ div = kdouble_div0_mod0(dvalue(tv_n), dvalue(tv_d), &mod);
+ tv_div = ktag_double(div);
+ tv_mod = ktag_double(mod);
+ break;
+ }
case K_TEINF:
if (ttiseinf(tv_n)) {
klispE_throw_simple(K, "non finite dividend");
@@ -1121,11 +1134,50 @@ void kdiv_mod(klisp_State *K, TValue *xparams, TValue ptree, TValue denv)
klispE_throw_simple(K, "non finite divisor");
return;
}
- default:
+ case K_TIINF:
+ if (ttisiinf(tv_n)) {
+ klispE_throw_simple(K, "non finite dividend");
+ return;
+ } else { /* if (ttiseinf(tv_d)) */
+ /* The semantics here are unclear, following the general
+ guideline of the report that says that if an infinity is
+ involved it should be understand as a limit. In that
+ case once the divisor is greater in magnitude than the
+ dividend the division stabilizes itself at q = 0; r = n
+ if both have the same sign, and q = 1; r = +infinity if
+ both have different sign (but in that case !(r < |d|)
+ !!) */
+ /* RATIONALE: if q were 0 we can't accomplish
+ q * d + r = n because q * d is undefined, if q isn't zero
+ then, either q*d + r is infinite or undefined so
+ there's no good q. on the other hand if we want
+ n - q*d = r & 0 <= r < d, r can't be infinite because it
+ would be equal to d, but q*d is infinite, so there's no
+ way out */
+ /* throw an exception, until this is resolved */
+ /* ASK John */
+ klispE_throw_simple(K, "non finite divisor");
+ return;
+ }
+ case K_TRWNPV: { /* no primary value */
+ /* ASK John: what happens with undefined & real with no primary values */
+ TValue n = ttisrwnpv(tv_n)? tv_n : tv_d;
+ if (kcurr_strict_arithp(K)) {
+ klispE_throw_simple_with_irritants(K, "operand has no primary "
+ "value", 1, n);
+ return;
+ } else {
+ tv_div = KRWNPV;
+ tv_mod = KRWNPV;
+ break;
+ }
+ }
+ default:
klispE_throw_simple(K, "unsupported type");
return;
}
+
TValue res;
if (flags & FDIV_DIV) {
if (flags & FDIV_MOD) { /* return both div and mod */
diff --git a/src/kreal.c b/src/kreal.c
@@ -656,3 +656,59 @@ void kdouble_print_string(klisp_State *K, TValue tv_double,
}
return;
}
+
+double kdouble_div_mod(double n, double d, double *res_mod)
+{
+ double div = floor(n / d);
+ double mod = n - div * d;
+
+ /* div, mod or div-and-mod */
+ /* 0 <= mod0 < |d| */
+ if (mod < 0.0) {
+ if (d < 0.0) {
+ mod -= d;
+ div += 1.0;
+ } else {
+ mod += d;
+ div -= 1.0;
+ }
+ }
+ *res_mod = mod;
+ return div;
+}
+
+double kdouble_div0_mod0(double n, double d, double *res_mod)
+{
+ double div = floor(n / d);
+ double mod = n - div * d;
+
+ /* div0, mod0 or div-and-mod0 */
+ /*
+ ** Adjust q and r so that:
+ ** -|d/2| <= mod0 < |d/2| which is the same as
+ ** dmin <= mod0 < dmax, where
+ ** dmin = -|d/2| and dmax = |d/2|
+ */
+ double dmin = -((d<0.0? -d : d) / 2.0);
+ double dmax = (d<0.0? -d : d) / 2.0;
+
+ if (mod < dmin) {
+ if (d < 0) {
+ mod -= d;
+ div += 1.0;
+ } else {
+ mod += d;
+ div -= 1.0;
+ }
+ } else if (mod >= dmax) {
+ if (d < 0) {
+ mod += d;
+ div += 1.0;
+ } else {
+ mod -= d;
+ div -= 1.0;
+ }
+ }
+ *res_mod = mod;
+ return div;
+}
diff --git a/src/kreal.h b/src/kreal.h
@@ -22,6 +22,9 @@ TValue kexact_to_inexact(klisp_State *K, TValue n);
TValue kinexact_to_exact(klisp_State *K, TValue n);
+double kdouble_div_mod(double n, double d, double *res_mod);
+double kdouble_div0_mod0(double n, double d, double *res_mod);
+
/*
** read/write interface
*/
