From fc14a0592d94b1dbb502afe298a67b9b544559af Mon Sep 17 00:00:00 2001
From: Kaz Kylheku <kaz@kylheku.com>
Date: Sat, 11 Jan 2014 23:09:22 -0800
Subject: * arith.c (rising_product): New static function. (n_choose_k,
 n_perm_k): New functions.

* arith.h (n_choose_k, n_perm_k): Declared.

* eval.c (eval_init): New functions interned.

* txr.1: Documented.
---
 ChangeLog | 11 +++++++++++
 arith.c   | 37 +++++++++++++++++++++++++++++++++++++
 arith.h   |  2 ++
 eval.c    |  2 ++
 txr.1     | 31 +++++++++++++++++++++++++++++++
 5 files changed, 83 insertions(+)

diff --git a/ChangeLog b/ChangeLog
index 40c1b7e2..14c4775a 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,14 @@
+2014-01-11  Kaz Kylheku  <kaz@kylheku.com>
+
+	* arith.c (rising_product): New static function.
+	(n_choose_k, n_perm_k): New functions.
+
+	* arith.h (n_choose_k, n_perm_k): Declared.
+
+	* eval.c (eval_init): New functions interned.
+
+	* txr.1: Documented.
+
 2014-01-11  Kaz Kylheku  <kaz@kylheku.com>
 
 	* arith.c (to_float): Print function name as ~a rather than ~s,
diff --git a/arith.c b/arith.c
index 445cebda..18a5fc37 100644
--- a/arith.c
+++ b/arith.c
@@ -1869,6 +1869,43 @@ val cum_norm_dist(val arg)
   }
 }
 
+static val rising_product(val m, val n)
+{
+  val acc;
+
+  if (lt(n, one))
+    return one;
+
+  if (ge(m, n))
+    return one;
+
+  if (lt(m, one))
+    m = one;
+
+  acc = m;
+
+  m = plus(m, one);
+
+  while (le(m, n)) {
+    acc = mul(acc, m);
+    m = plus(m, one);
+  }
+
+  return acc;
+}
+
+val n_choose_k(val n, val k)
+{
+  val top = rising_product(plus(minus(n, k), one), n);
+  val bottom = rising_product(one, k);
+  return trunc(top, bottom);
+}
+
+val n_perm_k(val n, val k)
+{
+  return rising_product(plus(minus(n, k), one), n);
+}
+
 void arith_init(void)
 {
   mp_init(&NUM_MAX_MP);
diff --git a/arith.h b/arith.h
index ef341707..a73bdde8 100644
--- a/arith.h
+++ b/arith.h
@@ -31,4 +31,6 @@ int highest_bit(int_ptr_t n);
 val normalize(val bignum);
 val in_int_ptr_range(val bignum);
 val cum_norm_dist(val x);
+val n_choose_k(val n, val k);
+val n_perm_k(val n, val k);
 void arith_init(void);
diff --git a/eval.c b/eval.c
index 12f3ebd0..03e8d34d 100644
--- a/eval.c
+++ b/eval.c
@@ -2353,6 +2353,8 @@ void eval_init(void)
   reg_fun(intern(lit("exp"), user_package), func_n1(expo));
   reg_fun(intern(lit("sqrt"), user_package), func_n1(sqroot));
   reg_fun(intern(lit("cum-norm-dist"), user_package), func_n1(cum_norm_dist));
+  reg_fun(intern(lit("n-choose-k"), user_package), func_n2(n_choose_k));
+  reg_fun(intern(lit("n-perm-k"), user_package), func_n2(n_perm_k));
   reg_fun(intern(lit("fixnump"), user_package), func_n1(fixnump));
   reg_fun(intern(lit("bignump"), user_package), func_n1(bignump));
   reg_fun(intern(lit("floatp"), user_package), func_n1(floatp));
diff --git a/txr.1 b/txr.1
index 20d871ef..a7812d15 100644
--- a/txr.1
+++ b/txr.1
@@ -8996,6 +8996,37 @@ The cum-norm-dist function calculates an approximation to the cumulative normal
 distribution function: the integral, of the normal distribution function, from
 negative infinity to the <argument>.
 
+.SS Functions n-choose-k and n-perm-k
+
+.TP
+Syntax:
+
+  (n-choose-k <n> <k>)
+  (n-perm-k <n> <k>)
+
+.TP
+Description:
+
+The n-choose-k function computes the binomial coefficient nCk which
+expresses the number of combinations of <k> items that can be chosen from
+a set of <n>, where combinations are subsets.
+
+The n-choose-k function computes nPk: the number of permutations of size <k>
+that can be drawn from a set of <n>, where permutations are sequences,
+whose order is significant.
+
+The calculations only make sense when <n> and <k> are nonnegative integers, and
+<k> does not exceed <n>. The behavior is not specified if these conditions
+are not met.
+
+.TP
+Description:
+
+The cum-norm-dist function calculates an approximation to the cumulative normal
+distribution function: the integral, of the normal distribution function, from
+negative infinity to the <argument>.
+
+
 .SS Functions fixnump, bignump, integerp, floatp, numberp
 
 .TP
-- 
cgit v1.2.3