4 files changed, 75 insertions, 0 deletions

diff --git a/arith.c b/arith.c
index 67fa75c2..cd3416ca 100644
--- a/arith.c
+++ b/arith.c
@@ -2102,6 +2102,59 @@ negop:
   uw_throw(error_s, lit("isqrt: negative operand"));
 }
 
+val square(val anum)
+{
+  val self = lit("square");
+
+  switch (type(anum)) {
+  case NUM:
+    {
+      cnum a = c_n(anum);
+#if HAVE_DOUBLE_INTPTR_T
+      double_intptr_t product = a * convert(double_intptr_t, a);
+      if (product < NUM_MIN || product > NUM_MAX)
+        return bignum_dbl_ipt(product);
+      return num_fast(product);
+#else
+      cnum ap = ABS(a);
+      if (2 * highest_bit(ap) < CNUM_BIT - 1) {
+        cnum product = a * a;
+        if (product >= NUM_MIN && product <= NUM_MAX)
+          return num_fast(product);
+        return bignum(product);
+      } else {
+        val n = make_bignum();
+        mp_err mpe;
+        mp_set_intptr(mp(n), a);
+        mpe = mp_sqr(mp(n), mp(n));
+        if (mpe != MP_OKAY)
+          do_mp_error(self, mpe);
+        return n;
+      }
+#endif
+    }
+  case BGNUM:
+    {
+      val n = make_bignum();
+      mp_err mpe = mp_sqr(mp(anum), mp(n));
+      if (mpe != MP_OKAY)
+        do_mp_error(self, mpe);
+      return n;
+    }
+  case FLNUM:
+    {
+      double a = c_flo(anum, self);
+      return flo(a * a);
+    }
+  case RNG:
+    return rcons(square(from(anum)), square(to(anum)));
+  default:
+    break;
+  }
+
+  uw_throwf(error_s, lit("square: invalid operand ~s"), anum, nao);
+}
+
 val gcd(val anum, val bnum)
 {
   val n;
diff --git a/eval.c b/eval.c
index 6ee34e7f..4985aeb9 100644
--- a/eval.c
+++ b/eval.c
@@ -6413,6 +6413,7 @@ void eval_init(void)
   reg_fun(intern(lit("expt"), user_package), func_n0v(exptv));
   reg_fun(intern(lit("exptmod"), user_package), func_n3(exptmod));
   reg_fun(intern(lit("isqrt"), user_package), func_n1(isqrt));
+  reg_fun(intern(lit("square"), user_package), func_n1(square));
   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_n2o(floordiv, 1));
diff --git a/lib.h b/lib.h
index 6570d904..d54ea7fa 100644
--- a/lib.h
+++ b/lib.h
@@ -721,6 +721,7 @@ val exptv(struct args *nlist);
 val exptmod(val base, val exp, val mod);
 val sqroot(val anum);
 val isqrt(val anum);
+val square(val anum);
 val gcd(val anum, val bnum);
 val gcdv(struct args *nlist);
 val lcm(val anum, val bnum);
diff --git a/txr.1 b/txr.1
index bf3c9353..f62cc5cf 100644
--- a/txr.1
+++ b/txr.1
@@ -35653,6 +35653,26 @@ and reduced to the
 least positive residue modulo
 .metn modulus .
 
+.coNP Function @ square
+.synb
+.mets (square << argument )
+.syne
+.desc
+The
+.code square
+function returns the product of
+.meta argument
+with itself. The following
+equivalence holds, except that
+.code x
+is evaluated only once in the the
+.code square
+expression:
+
+.cblk
+  (square x)  <-->  (* x x)
+.cble
+
 .coNP Function @ cum-norm-dist
 .synb
 .mets (cum-norm-dist << argument )