* arith.c (highest_bit): Linkage changed to static.

(abso, isqrt): New functions.
(isqrt_fixnum): New static function.

* eval.c (eval_init): Registered abs, sqrt and numberp instrinsics.

* lib.c (numberp): New function.

* lib.h (numberp, abso, isqrt): Declared.

* mpi-patches/series: New patch added.

* mpi-patches/faster-square-root: New patch added.

* txr.1: Documentation stubs for new functions.

8 files changed, 269 insertions, 3 deletions

diff --git a/ChangeLog b/ChangeLog
index 9fa763cb..0637746f 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,23 @@
 2011-12-13  Kaz Kylheku  <kaz@kylheku.com>
 
+	* arith.c (highest_bit): Linkage changed to static.
+	(abso, isqrt): New functions.
+	(isqrt_fixnum): New static function.
+
+	* eval.c (eval_init): Registered abs, sqrt and numberp instrinsics.
+
+	* lib.c (numberp): New function.
+
+	* lib.h (numberp, abso, isqrt): Declared.
+
+	* mpi-patches/series: New patch added.
+
+	* mpi-patches/faster-square-root: New patch added.
+
+	* txr.1: Documentation stubs for new functions.
+
+2011-12-13  Kaz Kylheku  <kaz@kylheku.com>
+
 	* arith.c (expt): Fix broken bignum x fixnum combination.
 
 2011-12-13  Kaz Kylheku  <kaz@kylheku.com>
diff --git a/arith.c b/arith.c
index 27034a2a..1ccab1f6 100644
--- a/arith.c
+++ b/arith.c
@@ -84,7 +84,7 @@ static val normalize(val bignum)
   }
 }
 
-int highest_bit(int_ptr_t n)
+static int highest_bit(int_ptr_t n)
 {
 #if SIZEOF_PTR == 8
   if (n & 0x7FFFFFFF00000000) {
@@ -414,6 +414,18 @@ val neg(val anum)
   }
 }
 
+val abso(val anum)
+{
+  if (bignump(anum)) {
+    val n = make_bignum();
+    mp_abs(mp(anum), mp(n));
+    return n;
+  } else {
+    cnum n = c_num(anum);
+    return num(n < 0 ? n : n);
+  }
+}
+
 val mul(val anum, val bnum)
 {
   int tag_a = tag(anum);
@@ -890,6 +902,36 @@ val expt(val anum, val bnum)
   abort();
 }
 
+static int_ptr_t isqrt_fixnum(int_ptr_t a)
+{
+  int_ptr_t mask = (int_ptr_t) 1 << (highest_bit(a) / 2);
+  int_ptr_t root = 0;
+
+  for (; mask != 0; mask >>= 1) {
+    int_ptr_t next_guess = root | mask;
+    if (next_guess * next_guess <= a)
+      root = next_guess;
+  }
+
+  return root;
+}
+
+val isqrt(val anum)
+{
+  if (fixnump(anum)) {
+    cnum a = c_num(anum);
+    if (a < 0)
+      uw_throw(error_s, lit("sqrt: negative operand"));
+    return num_fast(isqrt_fixnum(c_num(anum)));
+  } else if (bignump(anum)) {
+    val n = make_bignum();
+    if (mp_sqrt(mp(anum), mp(n)) != MP_OKAY)
+      uw_throw(error_s, lit("sqrt: negative operand"));
+    return normalize(n);
+  }
+  uw_throwf(error_s, lit("sqrt: invalid operand ~s"), anum, nao);
+}
+
 void arith_init(void)
 {
   mp_init(&NUM_MAX_MP);
diff --git a/eval.c b/eval.c
index 9c6cdf87..30a1beed 100644
--- a/eval.c
+++ b/eval.c
@@ -1156,11 +1156,14 @@ void eval_init(void)
   reg_fun(intern(lit("+"), user_package), func_n0v(plusv));
   reg_fun(intern(lit("-"), user_package), func_n1v(minusv));
   reg_fun(intern(lit("*"), user_package), func_n0v(mulv));
+  reg_fun(intern(lit("abs"), user_package), func_n1(abso));
   reg_fun(intern(lit("trunc"), user_package), func_n2(trunc));
   reg_fun(intern(lit("mod"), user_package), func_n2(mod));
   reg_fun(intern(lit("expt"), user_package), func_n0v(exptv));
+  reg_fun(intern(lit("sqrt"), user_package), func_n1(isqrt));
   reg_fun(intern(lit("fixnump"), user_package), func_n1(fixnump));
   reg_fun(intern(lit("bignump"), user_package), func_n1(bignump));
+  reg_fun(intern(lit("numberp"), user_package), func_n1(numberp));
 
   reg_fun(intern(lit("zerop"), user_package), func_n1(zerop));
   reg_fun(intern(lit(">"), user_package), func_n1v(gtv));
diff --git a/lib.c b/lib.c
index 21cfb68c..0062ac94 100644
--- a/lib.c
+++ b/lib.c
@@ -832,6 +832,20 @@ val bignump(val num)
   return (is_ptr(num) && type(num) == BGNUM) ? t : nil;
 }
 
+val numberp(val num)
+{
+  switch (tag(num)) {
+  case TAG_NUM:
+    return t;
+  case TAG_PTR:
+    if (num->t.type == BGNUM)
+      return t;
+    /* fallthrough */
+  default:
+    return nil;
+  }
+}
+
 val plusv(val nlist)
 {
   if (!nlist)
diff --git a/lib.h b/lib.h
index 4787e240..5874d73e 100644
--- a/lib.h
+++ b/lib.h
@@ -363,11 +363,13 @@ val num(cnum val);
 cnum c_num(val num);
 val fixnump(val num);
 val bignump(val num);
+val numberp(val num);
 val plus(val anum, val bnum);
 val plusv(val nlist);
 val minus(val anum, val bnum);
 val minusv(val minuend, val nlist);
 val neg(val num);
+val abso(val num);
 val mul(val anum, val bnum);
 val mulv(val nlist);
 val trunc(val anum, val bnum);
@@ -387,6 +389,7 @@ val maxv(val first, val rest);
 val minv(val first, val rest);
 val expt(val base, val exp);
 val exptv(val nlist);
+val isqrt(val anum);
 val string_own(wchar_t *str);
 val string(const wchar_t *str);
 val string_utf8(const char *str);
diff --git a/mpi-patches/faster-square-root b/mpi-patches/faster-square-root
new file mode 100644
index 00000000..a19bab03
--- /dev/null
+++ b/mpi-patches/faster-square-root
@@ -0,0 +1,183 @@
+Index: mpi-1.8.6/mpi.c
+===================================================================
+--- mpi-1.8.6.orig/mpi.c	2011-12-13 15:18:37.000000000 -0800
++++ mpi-1.8.6/mpi.c	2011-12-13 17:21:59.000000000 -0800
+@@ -158,6 +158,9 @@
+ mp_err   s_mp_grow(mp_int *mp, mp_size min);   /* increase allocated size */
+ mp_err   s_mp_pad(mp_int *mp, mp_size min);    /* left pad with zeroes    */
+ 
++int      s_highest_bit_mp(mp_int *a);
++mp_err   s_mp_set_bit(mp_int *a, int bit);
++
+ void     s_mp_clamp(mp_int *mp);               /* clip leading zeroes     */
+ 
+ void     s_mp_exch(mp_int *a, mp_int *b);      /* swap a and b in place   */
+@@ -1535,77 +1538,55 @@
+ 
+ /* {{{ mp_sqrt(a, b) */
+ 
+-/*
+-  mp_sqrt(a, b)
+-
+-  Compute the integer square root of a, and store the result in b.
+-  Uses an integer-arithmetic version of Newton's iterative linear
+-  approximation technique to determine this value; the result has the
+-  following two properties:
+-
+-     b^2 <= a
+-     (b+1)^2 >= a
+-
+-  It is a range error to pass a negative value.
+- */
+ mp_err mp_sqrt(mp_int *a, mp_int *b)
+ {
+-  mp_int   x, t;
+-  mp_err   res;
++  int mask_shift;
++  mp_int root, guess, *proot = &root, *pguess = &guess;
++  mp_int guess_sqr;
++  mp_err err = MP_MEM;
+ 
+   ARGCHK(a != NULL && b != NULL, MP_BADARG);
+ 
+-  /* Cannot take square root of a negative value */
+-  if(SIGN(a) == MP_NEG)
++  if (mp_cmp_z(b) == MP_LT)
+     return MP_RANGE;
+ 
+-  /* Special cases for zero and one, trivial     */
+-  if(mp_cmp_d(a, 0) == MP_EQ || mp_cmp_d(a, 1) == MP_EQ) 
+-    return mp_copy(a, b);
+-    
+-  /* Initialize the temporaries we'll use below  */
+-  if((res = mp_init_size(&t, USED(a))) != MP_OKAY)
+-    return res;
+-
+-  /* Compute an initial guess for the iteration as a itself */
+-  if((res = mp_init_copy(&x, a)) != MP_OKAY)
+-    goto X;
+-
+-  for(;;) {
+-    /* t = (x * x) - a */
+-    mp_copy(&x, &t);      /* can't fail, t is big enough for original x */
+-    if((res = mp_sqr(&t, &t)) != MP_OKAY ||
+-       (res = mp_sub(&t, a, &t)) != MP_OKAY)
+-      goto CLEANUP;
+-
+-    /* t = t / 2x       */
+-    s_mp_mul_2(&x);
+-    if((res = mp_div(&t, &x, &t, NULL)) != MP_OKAY)
+-      goto CLEANUP;
+-    s_mp_div_2(&x);
+-
+-    /* Terminate the loop, if the quotient is zero */
+-    if(mp_cmp_z(&t) == MP_EQ)
+-      break;
+-
+-    /* x = x - t       */
+-    if((res = mp_sub(&x, &t, &x)) != MP_OKAY)
+-      goto CLEANUP;
+-
++  if ((err = mp_init(&root)))
++    goto out;
++  if ((err = mp_init(&guess)))
++    goto cleanup_root;
++  if ((err = mp_init(&guess_sqr)))
++    goto cleanup_guess;
++
++  for (mask_shift = s_highest_bit_mp(a) / 2; mask_shift > 0; mask_shift--) {
++    mp_int *temp;
++
++    if ((err = mp_copy(proot, pguess)))
++	goto cleanup;
++    if ((err = s_mp_set_bit(pguess, mask_shift)))
++	goto cleanup;
++    if ((err = mp_copy(pguess, &guess_sqr)))
++	goto cleanup;
++    if ((err = s_mp_sqr(&guess_sqr)))
++	goto cleanup;
++
++    if (s_mp_cmp(&guess_sqr, a) <= 0) {
++	temp = proot;
++	proot = pguess;
++	pguess = temp;
++    }
+   }
+ 
+-  /* Copy result to output parameter */
+-  mp_sub_d(&x, 1, &x);
+-  s_mp_exch(&x, b);
+-
+- CLEANUP:
+-  mp_clear(&x);
+- X:
+-  mp_clear(&t); 
+-
+-  return res;
+-
+-} /* end mp_sqrt() */
++  err = mp_copy(proot, b);
++
++cleanup:
++  mp_clear(&guess_sqr);
++cleanup_guess:
++  mp_clear(&guess);
++cleanup_root:
++  mp_clear(&root);
++out:
++  return err;
++}
+ 
+ /* }}} */
+ 
+@@ -2541,21 +2522,9 @@
+ 
+ int    mp_count_bits(mp_int *mp)
+ {
+-  int      len;
+-  mp_digit d;
+-
+   ARGCHK(mp != NULL, MP_BADARG);
+ 
+-  len = DIGIT_BIT * (USED(mp) - 1);
+-  d = DIGIT(mp, USED(mp) - 1);
+-
+-  while(d != 0) {
+-    ++len;
+-    d >>= 1;
+-  }
+-
+-  return len;
+-  
++  return s_highest_bit_mp(mp);
+ } /* end mp_count_bits() */
+ 
+ /* }}} */
+@@ -3119,6 +3088,27 @@
+   abort();
+ }
+ 
++int s_highest_bit_mp(mp_int *a)
++{
++  int nd1 = USED(a) - 1;
++  return s_highest_bit(DIGIT(a, nd1)) + nd1 * MP_DIGIT_BIT;
++}
++
++mp_err s_mp_set_bit(mp_int *a, int bit)
++{
++  int nd = (bit + MP_DIGIT_BIT) / MP_DIGIT_BIT;
++  int nbit = bit - (nd - 1) * MP_DIGIT_BIT;
++  mp_err res;
++
++  if (nd == 0)
++    return MP_OKAY;
++
++  if ((res = s_mp_pad(a, nd)) != MP_OKAY)
++    return res;
++
++  DIGIT(a, nd - 1) |= ((mp_digit) 1 << nbit);
++  return MP_OKAY;
++}
+ 
+ /* {{{ s_mp_exch(a, b) */
+ 
diff --git a/mpi-patches/series b/mpi-patches/series
index 311d7fe5..0181c920 100644
--- a/mpi-patches/series
+++ b/mpi-patches/series
@@ -11,3 +11,4 @@ mpi-set-double-intptr
 fix-bad-shifts
 bit-search-optimizations
 shrink-mpi-int
+faster-square-root
diff --git a/txr.1 b/txr.1
index 8049fe55..a5ebdc06 100644
--- a/txr.1
+++ b/txr.1
@@ -4809,9 +4809,11 @@ The following are Lisp functions and variables built-in to TXR.
 
 .SS Functions eq, eql and equal
 
-.SS Arithmetic functions +, -, *, trunc, mod, expt
+.SS Arithmetic functions +, -, *, trunc, mod, expt, sqrt
 
-.SS Functions fixnump, bignump
+.SS Arithmetic function abs
+
+.SS Functions fixnump, bignump, numberp
 
 .SS Function zerop