floor and ceil do division, with optional second arg.

Also, with one argument, these functions handle
ranges.

* arith.c (floordiv, ceildiv): New functions.
(floorf, ceili): Handle ranges.

* eval.c (eval_init): Register floor and ceil
to new functions.

* lib.h (floordiv, ceildiv): Declared.

* txr.1: Documentation updated.

4 files changed, 228 insertions, 20 deletions

diff --git a/arith.c b/arith.c
index e1d105e8..3bded75b 100644
--- a/arith.c
+++ b/arith.c
@@ -1076,6 +1076,161 @@ val trunc_rem(val anum, val bnum)
   return list(quot, rem, nao);
 }
 
+val floordiv(val anum, val bnum)
+{
+  if (missingp(bnum))
+    return floorf(anum);
+tail:
+  switch (TAG_PAIR(tag(anum), tag(bnum))) {
+  case TAG_PAIR(TAG_NUM, TAG_NUM):
+    {
+      cnum a = c_num(anum);
+      cnum b = c_num(bnum);
+      cnum ap = ABS(a);
+      cnum bp = ABS(b);
+      int neg = ((a < 0 && b > 0) || (a > 0 && b < 0));
+
+      if (b == 0)
+        goto divzero;
+
+      {
+        cnum quot = ap / bp;
+        if (neg) {
+          if (quot * bp != ap)
+            return num(-quot - 1);
+          return num(-quot);
+        }
+        return num(quot);
+      }
+    }
+  case TAG_PAIR(TAG_NUM, TAG_PTR):
+    switch (type(bnum)) {
+    case BGNUM:
+      {
+        cnum a = c_num(anum);
+        if (a == 0)
+          return zero;
+        if (a < 0 && !ISNEG(mp(bnum)))
+          return negone;
+        if (a > 0 && ISNEG(mp(bnum)))
+          return negone;
+        return zero;
+      }
+    case FLNUM:
+      {
+        double x = c_num(anum), y = c_flo(bnum);
+        if (y == 0.0)
+          goto divzero;
+        else
+          return flo((x - dmod(x, y))/y);
+      }
+    default:
+      break;
+    }
+    break;
+  case TAG_PAIR(TAG_PTR, TAG_NUM):
+    switch (type(anum)) {
+    case BGNUM:
+      {
+        val n;
+        if (bnum == one)
+          return anum;
+        n = make_bignum();
+        if (sizeof (int_ptr_t) <= sizeof (mp_digit)) {
+          cnum b = c_num(bnum);
+          cnum bp = ABS(b);
+          mp_digit rem;
+          if (mp_div_d(mp(anum), bp, mp(n), &rem) != MP_OKAY)
+            goto divzero;
+          if (b < 0)
+            mp_neg(mp(n), mp(n));
+          if (rem && ((ISNEG(mp(anum)) && b > 0) ||
+                      (!ISNEG(mp(anum)) && b < 0)))
+            mp_sub_d(mp(n), 1, mp(n));
+        } else {
+          int err;
+          cnum b = c_num(bnum);
+          mp_int tmp, rem;
+          mp_init(&tmp);
+          mp_init(&rem);
+          mp_set_intptr(&tmp, b);
+          err = mp_div(mp(anum), &tmp, mp(n), 0);
+          mp_clear(&tmp);
+          if (err != MP_OKAY) {
+            mp_clear(&rem);
+            goto divzero;
+          }
+          if (mp_cmp_z(&rem) != MP_EQ &&
+              ((ISNEG(mp(anum)) && b > 0) ||
+               (!ISNEG(mp(anum)) && b < 0)))
+            mp_sub_d(mp(n), 1, mp(n));
+          mp_clear(&rem);
+        }
+        return normalize(n);
+      }
+    case FLNUM:
+      {
+        double x = c_flo(anum), y = c_num(bnum);
+        if (y == 0.0)
+          goto divzero;
+        else
+          return flo((x - dmod(x, y))/y);
+      }
+    case RNG:
+      return rcons(floordiv(from(anum), bnum), floordiv(to(anum), bnum));
+    default:
+      break;
+    }
+    break;
+  case TAG_PAIR(TAG_PTR, TAG_PTR):
+    switch (TYPE_PAIR(type(anum), type (bnum))) {
+    case TYPE_PAIR(BGNUM, BGNUM):
+      {
+        val n = make_bignum();
+        mp_int rem;
+        mp_init(&rem);
+        if (mp_div(mp(anum), mp(bnum), mp(n), &rem) != MP_OKAY) {
+          mp_clear(&rem);
+          goto divzero;
+        }
+        if (mp_cmp_z(&rem) != MP_EQ &&
+              ((ISNEG(mp(anum)) && !ISNEG(mp(bnum))) ||
+               (!ISNEG(mp(anum)) && ISNEG(mp(bnum)))))
+            mp_sub_d(mp(n), 1, mp(n));
+        mp_clear(&rem);
+        return normalize(n);
+      }
+    case TYPE_PAIR(FLNUM, FLNUM):
+      {
+        double x = c_flo(anum), y = c_flo(bnum);
+        if (y == 0.0)
+          goto divzero;
+        else
+          return flo((x - dmod(x, y))/y);
+      }
+    case TYPE_PAIR(BGNUM, FLNUM):
+      anum = flo_int(anum);
+      goto tail;
+    case TYPE_PAIR(FLNUM, BGNUM):
+      bnum = flo_int(bnum);
+      goto tail;
+    case TYPE_PAIR(RNG, BGNUM):
+    case TYPE_PAIR(RNG, FLNUM):
+      return rcons(floordiv(from(anum), bnum), floordiv(to(anum), bnum));
+    }
+  }
+  uw_throwf(error_s, lit("floor: invalid operands ~s ~s"), anum, bnum, nao);
+divzero:
+  uw_throw(numeric_error_s, lit("floor: division by zero"));
+}
+
+val ceildiv(val anum, val bnum)
+{
+  if (missingp(bnum))
+    return ceili(anum);
+  return neg(floordiv(neg(anum), bnum));
+}
+
 val wrap_star(val start, val end, val num)
 {
   val modulus = minus(end, start);
@@ -1688,16 +1843,34 @@ val lcm(val anum, val bnum)
 
 val floorf(val num)
 {
-  if (integerp(num))
+  switch (type(num)) {
+  case NUM:
+  case BGNUM:
     return num;
-  return flo(floor(c_flo(to_float(lit("floor"), num))));
+  case FLNUM:
+    return flo(floor(c_flo(num)));
+  case RNG:
+    return rcons(floorf(from(num)), floorf(to(num)));
+  default:
+    break;
+  }
+  uw_throwf(error_s, lit("floor: invalid operand ~s"), num);
 }
 
 val ceili(val num)
 {
-  if (integerp(num))
+  switch (type(num)) {
+  case NUM:
+  case BGNUM:
     return num;
-  return flo(ceil(c_flo(to_float(lit("ceil"), num))));
+  case FLNUM:
+    return flo(ceil(c_flo(num)));
+  case RNG:
+    return rcons(ceili(from(num)), ceili(to(num)));
+  default:
+    break;
+  }
+  uw_throwf(error_s, lit("ceil: invalid operand ~s"), num);
 }
 
 val sine(val num)
diff --git a/eval.c b/eval.c
index 418d79ff..41dd8e24 100644
--- a/eval.c
+++ b/eval.c
@@ -5702,8 +5702,8 @@ void eval_init(void)
   reg_fun(intern(lit("isqrt"), user_package), func_n1(isqrt));
   reg_fun(intern(lit("gcd"), user_package), func_n0v(gcdv));
   reg_fun(intern(lit("lcm"), user_package), func_n0v(lcmv));
-  reg_fun(intern(lit("floor"), user_package), func_n1(floorf));
-  reg_fun(intern(lit("ceil"), user_package), func_n1(ceili));
+  reg_fun(intern(lit("floor"), user_package), func_n2o(floordiv, 1));
+  reg_fun(intern(lit("ceil"), user_package), func_n2o(ceildiv, 1));
   reg_fun(intern(lit("sin"), user_package), func_n1(sine));
   reg_fun(intern(lit("cos"), user_package), func_n1(cosi));
   reg_fun(intern(lit("tan"), user_package), func_n1(tang));
diff --git a/lib.h b/lib.h
index 2e374e7c..6162c84f 100644
--- a/lib.h
+++ b/lib.h
@@ -679,7 +679,9 @@ val gcdv(struct args *nlist);
 val lcm(val anum, val bnum);
 val lcmv(struct args *nlist);
 val floorf(val);
+val floordiv(val, val);
 val ceili(val);
+val ceildiv(val anum, val bnum);
 val sine(val);
 val cosi(val);
 val tang(val);
diff --git a/txr.1 b/txr.1
index 6daca213..43cf2bf4 100644
--- a/txr.1
+++ b/txr.1
@@ -31690,27 +31690,60 @@ returned: a positive number of the same type with exactly the same magnitude.
 
 .coNP Functions @ floor and @ ceil
 .synb
-.mets (floor << number )
-.mets (ceil << number )
+.mets (floor < dividend <> [ divisor ])
+.mets (ceil < dividend <> [ divisor ])
 .syne
 .desc
 The
 .code floor
+and
+.code ceiling
+functions perform division of the
+.meta dividend
+by the
+.metn divisor ,
+returning an integer quotient.
+
+If the
+.meta divisor
+is omitted, it defaults to 1.
+
+If both inputs are integers,
+the result is of type integer.
+
+If all inputs are numbers and at least one of them is
+floating-point, the others are converted to floating-point
+and the result is floating-point.
+
+The
+.code dividend
+input may be a range. In this situation, the operation is
+recursively distributed over the
+.code from
+and
+.code to
+fields of the range, individually matched against the
+.metn divisor ,
+and the result is a range composed of these two individual
+quotients.
+
+When the quotient is a scalar value,
+.code floor
 function returns the highest integer which does not exceed
-the value of
-.metn number .
-The ceiling function returns the lowest integer which
+the value of the quotient. That is to say, the division is
+truncated to an integer value toward negative infinity.
+The
+.code ceil
+function the lowest integer which is not below the value
+of the quotient.
 does not exceed the value of
-.metn number .
-
-If
-.meta number
-an integer, it is simply returned.
+.metn dividend .
+That is to say, the division is truncated to an integer
+value toward positive infinity.
 
-If the argument is a float, then the value returned is a float.
-For instance
-.code "(floor 1.1)"
-returns 1.0 rather than 1.
+Note that for large floating point values, due to the limited
+precision, the integer value corresponding to the mathematical
+floor or ceiling may not be available.
 
 .coNP Functions @, sin @, cos @, tan @, asin @, acos @ atan and @ atan2
 .synb