From a12728d5a310a2964276823ffab0cb42f224c698 Mon Sep 17 00:00:00 2001
From: John Benediktsson <mrjbq7@gmail.com>
Date: Mon, 8 Apr 2013 15:03:15 -0700
Subject: [PATCH] math.factorials: various factorial words.

---
 extra/math/factorials/authors.txt             |  1 +
 extra/math/factorials/factorials-tests.factor | 36 +++++++++++++
 extra/math/factorials/factorials.factor       | 53 +++++++++++++++++++
 3 files changed, 90 insertions(+)
 create mode 100644 extra/math/factorials/authors.txt
 create mode 100644 extra/math/factorials/factorials-tests.factor
 create mode 100644 extra/math/factorials/factorials.factor

diff --git a/extra/math/factorials/authors.txt b/extra/math/factorials/authors.txt
new file mode 100644
index 0000000000..e091bb8164
--- /dev/null
+++ b/extra/math/factorials/authors.txt
@@ -0,0 +1 @@
+John Benediktsson
diff --git a/extra/math/factorials/factorials-tests.factor b/extra/math/factorials/factorials-tests.factor
new file mode 100644
index 0000000000..cd110d515a
--- /dev/null
+++ b/extra/math/factorials/factorials-tests.factor
@@ -0,0 +1,36 @@
+USING: kernel math.functions tools.test ;
+IN: math.factorials
+
+[ 1 ] [ -1 factorial ] unit-test ! not necessarily correct
+[ 1 ] [ 0 factorial ] unit-test
+[ 1 ] [ 1 factorial ] unit-test
+[ 3628800 ] [ 10 factorial ] unit-test
+
+{ 1 } [ 10 10 factorial/ ] unit-test
+{ 720 } [ 10 7 factorial/ ] unit-test
+{ 604800 } [ 10 3 factorial/ ] unit-test
+{ 3628800 } [ 10 0 factorial/ ] unit-test
+
+{ 17160 } [ 10 4 rising-factorial ] unit-test
+{ 1/57120 } [ 10 -4 rising-factorial ] unit-test
+{ 10 } [ 10 1 rising-factorial ] unit-test
+{ 0 } [ 10 0 rising-factorial ] unit-test
+
+{ 5040 } [ 10 4 falling-factorial ] unit-test
+{ 1/24024 } [ 10 -4 falling-factorial ] unit-test
+{ 10 } [ 10 1 falling-factorial ] unit-test
+{ 0 } [ 10 0 falling-factorial ] unit-test
+
+{ 7301694400 } [ 100 5 3 factorial-power ] unit-test
+{ 5814000000 } [ 100 5 5 factorial-power ] unit-test
+{ 4549262400 } [ 100 5 7 factorial-power ] unit-test
+{ 384000000 } [ 100 5 20 factorial-power ] unit-test
+{ 384000000 } [ 100 5 20 factorial-power ] unit-test
+{ 44262400 } [ 100 5 24 factorial-power ] unit-test
+{ 0 } [ 100 5 25 factorial-power ] unit-test
+{ 4760 } [ 20 3 3 factorial-power ] unit-test
+{ 1/17342 } [ 20 -3 3 factorial-power ] unit-test
+{ 1/2618 } [ 20 -3 -3 factorial-power ] unit-test
+{ 11960 } [ 20 3 -3 factorial-power ] unit-test
+{ t } [ 20 3 [ 1 factorial-power ] [ falling-factorial ] 2bi = ] unit-test
+{ t } [ 20 3 [ 0 factorial-power ] [ ^ ] 2bi = ] unit-test
diff --git a/extra/math/factorials/factorials.factor b/extra/math/factorials/factorials.factor
new file mode 100644
index 0000000000..07664fc89b
--- /dev/null
+++ b/extra/math/factorials/factorials.factor
@@ -0,0 +1,53 @@
+! Copyright (C) 2013 John Benediktsson
+! See http://factorcode.org/license.txt for BSD license
+
+USING: combinators combinators.short-circuit kernel locals math
+math.functions math.ranges memoize sequences ;
+
+IN: math.factorials
+
+MEMO: factorial ( n -- n! )
+    dup 1 > [ [1,b] product ] [ drop 1 ] if ;
+
+:: factorial/ ( n k -- n!/k! )
+    { [ k 0 < ] [ n 0 < ] [ k n > ] } 0||
+    [ 0 ] [ k n (a,b] product ] if ;
+
+: rising-factorial ( x n -- x(n) )
+    {
+        { 1 [ ] }
+        { 0 [ drop 0 ] }
+        [
+            dup 0 < [ neg [ + ] keep t ] [ f ] if
+            [ dupd + [a,b) product ] dip
+            [ recip ] when
+        ]
+    } case ;
+
+ALIAS: pochhammer rising-factorial
+
+: falling-factorial ( x n -- (x)n )
+    {
+        { 1 [ ] }
+        { 0 [ drop 0 ] }
+        [
+            dup 0 < [ neg [ + ] keep t ] [ f ] if
+            [ dupd - swap (a,b] product ] dip
+            [ recip ] when
+        ]
+    } case ;
+
+: factorial-power ( x n h -- (x)n(h) )
+    {
+        { 1 [ falling-factorial ] }
+        { 0 [ ^ ] }
+        [
+            over 0 < [
+                [ [ nip + ] [ swap neg * + ] 3bi ] keep
+                <range> product recip
+            ] [
+                neg [ [ dupd 1 - ] [ * ] bi* + ] keep
+                <range> product
+            ] if
+        ]
+    } case ;