From 55525a8ba48943e279adde19cd79c92b607b6ab8 Mon Sep 17 00:00:00 2001
From: John Benediktsson <mrjbq7@gmail.com>
Date: Fri, 1 Jun 2012 14:07:39 -0700
Subject: [PATCH] math.combinatorics: 100%+ faster combinations.

---
 .../combinatorics/combinatorics-tests.factor  | 20 +-----
 basis/math/combinatorics/combinatorics.factor | 72 +++++++------------
 2 files changed, 29 insertions(+), 63 deletions(-)

diff --git a/basis/math/combinatorics/combinatorics-tests.factor b/basis/math/combinatorics/combinatorics-tests.factor
index 4e61f27f0d..5ae2ceffc0 100644
--- a/basis/math/combinatorics/combinatorics-tests.factor
+++ b/basis/math/combinatorics/combinatorics-tests.factor
@@ -52,24 +52,8 @@ IN: math.combinatorics.tests
 [ "ABC" "ACB" "BAC" "BCA" "CAB" "CBA" "ABC" ]
 [ "ABC" 6 [ dup dup clone-like next-permutation ] times ] unit-test
 
-[ 2598960 ] [ 52 iota 5 <combo> choose ] unit-test
-
-[ 6 3 13 6 ] [ 7 4 28 next-values ] unit-test
-[ 5 2 3 5 ] [ 6 3 13 next-values ] unit-test
-[ 3 1 0 3 ] [ 5 2 3 next-values ] unit-test
-[ 0 0 0 0 ] [ 3 1 0 next-values ] unit-test
-
-[ 9 ] [ 0 5 iota 3 <combo> dual-index ] unit-test
-[ 0 ] [ 9 5 iota 3 <combo> dual-index ] unit-test
-[ 179 ] [ 72 10 iota 5 <combo> dual-index ] unit-test
-
-[ { 5 3 2 1 } ] [ 7 iota 4 <combo> 8 combinadic ] unit-test
-[ { 4 3 2 1 0 } ] [ 10 iota 5 <combo> 0 combinadic ] unit-test
-[ { 8 6 3 1 0 } ] [ 10 iota 5 <combo> 72 combinadic ] unit-test
-[ { 9 8 7 6 5 } ] [ 10 iota 5 <combo> 251 combinadic ] unit-test
-
-[ { 0 1 2 } ] [ 0 5 iota 3 <combo> combination-indices ] unit-test
-[ { 2 3 4 } ] [ 9 5 iota 3 <combo> combination-indices ] unit-test
+[ { 0 1 2 } ] [ 0 3 5 combination-indices ] unit-test
+[ { 2 3 4 } ] [ 9 3 5 combination-indices ] unit-test
 
 [ { "a" "b" "c" } ] [ 0 { "a" "b" "c" "d" "e" } 3 combination ] unit-test
 [ { "c" "d" "e" } ] [ 9 { "a" "b" "c" "d" "e" } 3 combination ] unit-test
diff --git a/basis/math/combinatorics/combinatorics.factor b/basis/math/combinatorics/combinatorics.factor
index fd0ad7ba7b..665629dce7 100644
--- a/basis/math/combinatorics/combinatorics.factor
+++ b/basis/math/combinatorics/combinatorics.factor
@@ -4,7 +4,7 @@
 USING: accessors arrays assocs binary-search classes.tuple fry
 kernel locals math math.order math.ranges namespaces sequences
 sequences.private sorting ;
-
+FROM: sequences => change-nth ;
 IN: math.combinatorics
 
 <PRIVATE
@@ -113,61 +113,43 @@ PRIVATE>
 
 <PRIVATE
 
-TUPLE: combo
-    { seq sequence }
-    { k integer } ;
+! "Algorithm 515: Generation of a Vector from the Lexicographical Index"
+! Buckles, B. P., and Lybanon, M. ACM
+! Transactions on Mathematical Software, Vol. 3, No. 2, June 1977.
 
-C: <combo> combo
+:: combination-indices ( x! p n -- seq )
+    x 1 + x!
+    p 0 <array> :> c 0 :> k! 0 :> r!
+    p 1 - [| i |
+        i [ 0 ] [ 1 - c nth ] if-zero i c set-nth
+        [ k x < ] [
+            i c [ 1 + ] change-nth
+            n i c nth - p i 1 + - nCk r!
+            k r + k!
+        ] do while k r - k!
+    ] each-integer
+    p 2 < [ 0 ] [ p 2 - c nth ] if
+    p 1 < [ drop ] [ x + k - p 1 - c set-nth ] if
+    c [ 1 - ] map! ;
 
-: choose ( combo -- nCk )
-    [ seq>> length ] [ k>> ] bi nCk ;
-
-: largest-value ( a b x -- v )
-    dup 0 = [
-        drop 1 - nip
-    ] [
-        [ iota ] 2dip '[ _ nCk _ >=< ] search nip
-    ] if ;
-
-:: next-values ( a b x -- a' b' x' v )
-    a b x largest-value dup :> v  ! a'
-    b 1 -                         ! b'
-    x v b nCk -                   ! x'
-    v ;                           ! v == a'
-
-: dual-index ( m combo -- m' )
-    choose 1 - swap - ;
-
-: initial-values ( combo m -- n k m )
-    [ [ seq>> length ] [ k>> ] bi ] dip ;
-
-: combinadic ( combo m -- combinadic )
-    initial-values [ over 0 > ] [ next-values ] produce
-    [ 3drop ] dip ;
-
-:: combination-indices ( m combo -- seq )
-    combo m combo dual-index combinadic
-    combo seq>> length 1 - swap [ - ] with map! ;
-
-: apply-combination ( m combo -- seq )
-    [ combination-indices ] keep seq>> nths ;
-
-: combinations-quot ( seq k quot -- seq quot )
-    [ <combo> [ choose iota ] keep ] dip
-    '[ _ apply-combination @ ] ; inline
+:: combinations-quot ( seq k quot -- seq quot )
+    seq length :> n
+    n k nCk iota [
+        k n combination-indices seq nths quot call
+    ] ; inline
 
 PRIVATE>
 
 : combination ( m seq k -- seq' )
-    <combo> apply-combination ;
+    swap [ length combination-indices ] [ nths ] bi ;
 
-TUPLE: combinations combo length ;
+TUPLE: combinations seq k length ;
 
 : <combinations> ( seq k -- combinations )
-    [ <combo> ] 2keep [ length ] [ nCk ] bi* combinations boa ;
+    2dup [ length ] [ nCk ] bi* combinations boa ;
 
 M: combinations length length>> ; inline
-M: combinations nth-unsafe combo>> apply-combination ;
+M: combinations nth-unsafe [ seq>> ] [ k>> ] bi combination ;
 M: combinations hashcode* tuple-hashcode ;
 
 INSTANCE: combinations immutable-sequence