VBAHaskellで順列の数え上げ

n 個の元から k-個を選んで得られる順列の総数は下の式で計算できる。

VBAの適当な配列から、このすべての並びを列挙する関数 nPk をVBAHaskellで書いてみる。
（CallMeKohei さんのつぶやきからネタをパクりましたm(_ _ )m）

関数 nPk に任意の配列と取り出す要素数（省略可）を与えると、取りうる全ての並びからなる配列を返す。それぞれの要素それ自体もまた配列なので、返り値はジャグ配列ということになる。5つの文字からなる配列 Array("a", "b", "c", "d", "e") で結果を確認すると次のようになる。

m = nPk(Array("a", "b", "c", "d", "e"))    ' 全体の順列
?sizeof(m)
 120 
printM unzip(m,2)
  a  b  c  d  e
  a  b  c  e  d
 ～～省略～～～
  e  d  c  a  b
  e  d  c  b  a

m = nPk(Array("a", "b", "c", "d", "e"), 2)    ' 2個だけ取り出す順列
?sizeof(m)
 20 
printM unzip(m,2)
  a  b
  a  c
 ～～省略～～～
  e  c
  e  d

与えた配列に重複する要素があった場合には取り除かなければならないが、それは呼び出し側の責任とした。
VBAHaskellを使った nPk 関数のソースコードは以下の通り。

    ' 順列の数え上げ
' ar :  任意の1次元配列  ,  k_val :  取り出す個数（省略時は配列全体）
Function nPk(ByRef ar As Variant, Optional ByRef k_val As Variant)
    If IsMissing(k_val) Then k_val = sizeof(ar)
    If sizeof(ar) < k_val Then k_val = sizeof(ar)
    If 1 < k_val Then
        Dim ret As Variant: ret = VBA.Array()
        Dim flg As Variant: flg = iota(LBound(ar), UBound(ar))
        Dim i As Long
        ' 1要素を選択　→　その要素と残りの要素に分割　→　残り要素に対して再帰　→　つなげる
        For i = LBound(ar) To UBound(ar) Step 1
            ret = catV(ret, mapF(p_cons(ar(i)), nPk(filterR(ar, mapF(p_notEqual(i), flg)), k_val - 1)))
        Next i
        swapVariant nPk, ret
    ElseIf 1 = k_val Then
        nPk = ar
    Else
        nPk = VBA.Array()
    End If
End Function

このプログラムにはライブラリとしてVBAHaskellが必要で、ここで使用している主な関数を表にまとめた。 qiita.com

関数	機能	配置モジュール	備考
mapF	配列の各要素に関数を適用	Haskell_1_Core
notEqual	述語 a != b	Haskell_2_stdFun	関数オブジェクト p_notEqual として使用
filterR	配列要素のフィルタリング	Haskell_4_vector
cons	配列先頭への要素追加	Haskell_4_vector	関数オブジェクト p_cons として使用
catV	配列の結合	Haskell_4_vector
sizeof	配列長の取得	Haskell_4_vector
iota	連続する整数列の生成	Haskell_4_vector
swapVariant	Variant変数のスワップ	Haskell_0_declare

これに基づいていちばん複雑な式

ret = catV(ret, mapF(p_cons(ar(i)), nPk(filterR(ar, mapF(p_notEqual(i), flg)), k_val - 1)))

の説明をする。その前の段階で flg = iota(LBound(ar), UBound(ar))としているので flg には配列のインデックスが入っている。ar = Array("a", "b", "c", "d", "e") だったら flg = Array(0, 1, 2, 3, 4) だ。ループインデックスiがカウントアップされたときの ar(i) と mapF(p_notEqual(i), flg) と filterR(ar, mapF(p_notEqual(i), flg)) の値は以下のように推移する。

i	ar(i)	mapF(p_notEqual(i), flg)	filterR(ar, mapF(p_notEqual(i), flg))
0	"a"	Array(0, 1, 1, 1, 1)	Array("b", "c", "d", "e")
1	"b"	Array(1, 0, 1, 1, 1)	Array("a", "c", "d", "e")
2	"c"	Array(1, 1, 0, 1, 1)	Array("a", "b", "d", "e")
3	"d"	Array(1, 1, 1, 0, 1)	Array("a", "b", "c", "e")
4	"e"	Array(1, 1, 1, 1, 0)	Array("a", "b", "c", "d")

だから再帰呼び出し nPk(filterR(ar, mapF(p_notEqual(i), flg)), k_val - 1) は先頭1要素を除いた残りの配列に対する再帰になっている。返ってくる値はジャグ配列となるがそれぞれの先頭に最初選んだ1要素を付け加える、というのをくりかえす。

配列の先頭に要素を付け加える関数は cons なので、その第1引数を ar(i) に束縛して配列にマップしてやればいい。式 p_cons(ar(i)) がそれである。
一番外側にある catV 関数は配列の結合で、先頭 "a" のグループ、先頭 "b" のグループ・・・先頭 "e" のグループを結合してひとつの配列にしている。

ただしこの関数は実用性がない。手元のマシンで要素数7のとき（順列の数 = 5040）かかった時間は約1秒、要素数8のとき（順列の数 = 40320）かかった時間は約9秒だ。