--- /dev/null
+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
--- /dev/null
+! 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 ;