import itertools
_pascal_triangle = {(0, 0): 1}
def choose(n, r):
"""Binomial coefficient nCr.
Implemented using a cache. O(1) on cache hit, up to O(n*r) on miss.
"""
if r < 0 or r > n:
return 0
key = (n, r)
res = _pascal_triangle.get(key)
if res is not None:
return res
res = choose(n - 1, r - 1) + choose(n - 1, r)
_pascal_triangle[key] = res
return res
def first_naive(i, n, r):
"""Find first element and index of first combination with that first element.
Returns a tuple of value and index.
Naive implementation using an O(n) loop, with one choose invocation per loop.
Example: first_naive(8, 5, 3) returns (1, 6) because the combination with
index 8 is [1, 3, 4] so it starts with 1, and because the first combination
that starts with 1 is [1, 2, 3] which has index 6.
"""
s1 = 0
for k in range(n):
s2 = s1 + choose(n - k - 1, r - 1)
if i < s2:
return k, s1
s1 = s2
def first_bisect(i, n, r):
"""Find first element and index of first combination with that first element.
Returns a tuple of value and index.
Efficient implementation using O(log n) bisection steps, with one choose
invocation in each step.
Example: first_bisect(8, 5, 3) returns (1, 6) because the combination with
index 8 is [1, 3, 4] so it starts with 1, and because the first combination
that starts with 1 is [1, 2, 3] which has index 6.
"""
nCr = choose(n, r)
k1 = r - 1
s1 = nCr
k2 = n
s2 = 0
while k2 - k1 > 1:
k3 = (k1 + k2) // 2
s3 = nCr - choose(k3, r)
if s3 <= i:
k2, s2 = k3, s3
else:
k1, s1 = k3, s3
return n - k2, s2
def combination(i, n, r):
"""Compute combination with a given index.
Equivalent to list(itertools.combinations(range(n), r))[i].
Each combination is represented as a tuple of ascending elements, and
combinations are ordered lexicograplically.
Args:
i: zero-based index of the combination
n: number of possible values, will be taken from range(n)
r: number of elements in result list
"""
if r == 0:
return []
k, ik = first_bisect(i, n, r)
return tuple([k] + [j + k + 1 for j in combination(i - ik, n - k - 1, r - 1)])
def show_all(n, r):
"""Print all combinations for a given n and r, in order."""
print("n={} r={}".format(n, r))
for i in range(choose(n, r)):
print(" {:8d}: {}".format(i, combination(i, n, r)))
print("")
def test_first():
"""Check that first)naive and first_bisect are equivalent."""
for n in range(2, 10):
for r in range(1, n):
for i in range(choose(n, r)):
want = first_naive(i, n, r)
got = first_bisect(i, n, r)
assert want == got, "i={} n={} r={}: {} != {}".format(i, n, r, want, got)
def test_combination():
"""Check that combination and itertools.combinations are equivalent."""
for n in range(2, 10):
for r in range(1, n):
got = [combination(i, n, r) for i in range(choose(n, r))]
want = list(itertools.combinations(range(n), r))
assert want == got, "n={} r={}: {} != {}".format(n, r, want, got)
show_all(5, 3)
test_first()
test_combination()
aW1wb3J0IGl0ZXJ0b29scwoKX3Bhc2NhbF90cmlhbmdsZSA9IHsoMCwgMCk6IDF9CgpkZWYgY2hvb3NlKG4sIHIpOgogICIiIkJpbm9taWFsIGNvZWZmaWNpZW50IG5Dci4KCiAgSW1wbGVtZW50ZWQgdXNpbmcgYSBjYWNoZS4gTygxKSBvbiBjYWNoZSBoaXQsIHVwIHRvIE8obipyKSBvbiBtaXNzLgogICIiIgogIGlmIHIgPCAwIG9yIHIgPiBuOgogICAgcmV0dXJuIDAKICBrZXkgPSAobiwgcikKICByZXMgPSBfcGFzY2FsX3RyaWFuZ2xlLmdldChrZXkpCiAgaWYgcmVzIGlzIG5vdCBOb25lOgogICAgcmV0dXJuIHJlcwogIHJlcyA9IGNob29zZShuIC0gMSwgciAtIDEpICsgY2hvb3NlKG4gLSAxLCByKQogIF9wYXNjYWxfdHJpYW5nbGVba2V5XSA9IHJlcwogIHJldHVybiByZXMKCmRlZiBmaXJzdF9uYWl2ZShpLCBuLCByKToKICAiIiJGaW5kIGZpcnN0IGVsZW1lbnQgYW5kIGluZGV4IG9mIGZpcnN0IGNvbWJpbmF0aW9uIHdpdGggdGhhdCBmaXJzdCBlbGVtZW50LgoKICBSZXR1cm5zIGEgdHVwbGUgb2YgdmFsdWUgYW5kIGluZGV4LgoKICBOYWl2ZSBpbXBsZW1lbnRhdGlvbiB1c2luZyBhbiBPKG4pIGxvb3AsIHdpdGggb25lIGNob29zZSBpbnZvY2F0aW9uIHBlciBsb29wLgoKICBFeGFtcGxlOiBmaXJzdF9uYWl2ZSg4LCA1LCAzKSByZXR1cm5zICgxLCA2KSBiZWNhdXNlIHRoZSBjb21iaW5hdGlvbiB3aXRoCiAgaW5kZXggOCBpcyBbMSwgMywgNF0gc28gaXQgc3RhcnRzIHdpdGggMSwgYW5kIGJlY2F1c2UgdGhlIGZpcnN0IGNvbWJpbmF0aW9uCiAgdGhhdCBzdGFydHMgd2l0aCAxIGlzIFsxLCAyLCAzXSB3aGljaCBoYXMgaW5kZXggNi4KICAiIiIKICBzMSA9IDAKICBmb3IgayBpbiByYW5nZShuKToKICAgIHMyID0gczEgKyBjaG9vc2UobiAtIGsgLSAxLCByIC0gMSkKICAgIGlmIGkgPCBzMjoKICAgICAgcmV0dXJuIGssIHMxCiAgICBzMSA9IHMyCgpkZWYgZmlyc3RfYmlzZWN0KGksIG4sIHIpOgogICIiIkZpbmQgZmlyc3QgZWxlbWVudCBhbmQgaW5kZXggb2YgZmlyc3QgY29tYmluYXRpb24gd2l0aCB0aGF0IGZpcnN0IGVsZW1lbnQuCgogIFJldHVybnMgYSB0dXBsZSBvZiB2YWx1ZSBhbmQgaW5kZXguCgogIEVmZmljaWVudCBpbXBsZW1lbnRhdGlvbiB1c2luZyBPKGxvZyBuKSBiaXNlY3Rpb24gc3RlcHMsIHdpdGggb25lIGNob29zZQogIGludm9jYXRpb24gaW4gZWFjaCBzdGVwLgoKICBFeGFtcGxlOiBmaXJzdF9iaXNlY3QoOCwgNSwgMykgcmV0dXJucyAoMSwgNikgYmVjYXVzZSB0aGUgY29tYmluYXRpb24gd2l0aAogIGluZGV4IDggaXMgWzEsIDMsIDRdIHNvIGl0IHN0YXJ0cyB3aXRoIDEsIGFuZCBiZWNhdXNlIHRoZSBmaXJzdCBjb21iaW5hdGlvbgogIHRoYXQgc3RhcnRzIHdpdGggMSBpcyBbMSwgMiwgM10gd2hpY2ggaGFzIGluZGV4IDYuCiAgIiIiCiAgbkNyID0gY2hvb3NlKG4sIHIpCiAgazEgPSByIC0gMQogIHMxID0gbkNyCiAgazIgPSBuCiAgczIgPSAwCiAgd2hpbGUgazIgLSBrMSA+IDE6CiAgICBrMyA9IChrMSArIGsyKSAvLyAyCiAgICBzMyA9IG5DciAtIGNob29zZShrMywgcikKICAgIGlmIHMzIDw9IGk6CiAgICAgIGsyLCBzMiA9IGszLCBzMwogICAgZWxzZToKICAgICAgazEsIHMxID0gazMsIHMzCiAgcmV0dXJuIG4gLSBrMiwgczIKCmRlZiBjb21iaW5hdGlvbihpLCBuLCByKToKICAiIiJDb21wdXRlIGNvbWJpbmF0aW9uIHdpdGggYSBnaXZlbiBpbmRleC4KCiAgRXF1aXZhbGVudCB0byBsaXN0KGl0ZXJ0b29scy5jb21iaW5hdGlvbnMocmFuZ2UobiksIHIpKVtpXS4KCiAgRWFjaCBjb21iaW5hdGlvbiBpcyByZXByZXNlbnRlZCBhcyBhIHR1cGxlIG9mIGFzY2VuZGluZyBlbGVtZW50cywgYW5kCiAgY29tYmluYXRpb25zIGFyZSBvcmRlcmVkIGxleGljb2dyYXBsaWNhbGx5LgoKICBBcmdzOgogICAgaTogemVyby1iYXNlZCBpbmRleCBvZiB0aGUgY29tYmluYXRpb24KICAgIG46IG51bWJlciBvZiBwb3NzaWJsZSB2YWx1ZXMsIHdpbGwgYmUgdGFrZW4gZnJvbSByYW5nZShuKQogICAgcjogbnVtYmVyIG9mIGVsZW1lbnRzIGluIHJlc3VsdCBsaXN0CiAgIiIiCiAgaWYgciA9PSAwOgogICAgcmV0dXJuIFtdCiAgaywgaWsgPSBmaXJzdF9iaXNlY3QoaSwgbiwgcikKICByZXR1cm4gdHVwbGUoW2tdICsgW2ogKyBrICsgMSBmb3IgaiBpbiBjb21iaW5hdGlvbihpIC0gaWssIG4gLSBrIC0gMSwgciAtIDEpXSkKCmRlZiBzaG93X2FsbChuLCByKToKICAiIiJQcmludCBhbGwgY29tYmluYXRpb25zIGZvciBhIGdpdmVuIG4gYW5kIHIsIGluIG9yZGVyLiIiIgogIHByaW50KCJuPXt9IHI9e30iLmZvcm1hdChuLCByKSkKICBmb3IgaSBpbiByYW5nZShjaG9vc2UobiwgcikpOgogICAgcHJpbnQoIiAgezo4ZH06IHt9Ii5mb3JtYXQoaSwgY29tYmluYXRpb24oaSwgbiwgcikpKQogIHByaW50KCIiKQoKZGVmIHRlc3RfZmlyc3QoKToKICAiIiJDaGVjayB0aGF0IGZpcnN0KW5haXZlIGFuZCBmaXJzdF9iaXNlY3QgYXJlIGVxdWl2YWxlbnQuIiIiCiAgZm9yIG4gaW4gcmFuZ2UoMiwgMTApOgogICAgZm9yIHIgaW4gcmFuZ2UoMSwgbik6CiAgICAgIGZvciBpIGluIHJhbmdlKGNob29zZShuLCByKSk6CiAgICAgICAgd2FudCA9IGZpcnN0X25haXZlKGksIG4sIHIpCiAgICAgICAgZ290ID0gZmlyc3RfYmlzZWN0KGksIG4sIHIpCiAgICAgICAgYXNzZXJ0IHdhbnQgPT0gZ290LCAiaT17fSBuPXt9IHI9e306IHt9ICE9IHt9Ii5mb3JtYXQoaSwgbiwgciwgd2FudCwgZ290KQoKZGVmIHRlc3RfY29tYmluYXRpb24oKToKICAiIiJDaGVjayB0aGF0IGNvbWJpbmF0aW9uIGFuZCBpdGVydG9vbHMuY29tYmluYXRpb25zIGFyZSBlcXVpdmFsZW50LiIiIgogIGZvciBuIGluIHJhbmdlKDIsIDEwKToKICAgIGZvciByIGluIHJhbmdlKDEsIG4pOgogICAgICBnb3QgPSBbY29tYmluYXRpb24oaSwgbiwgcikgZm9yIGkgaW4gcmFuZ2UoY2hvb3NlKG4sIHIpKV0KICAgICAgd2FudCA9IGxpc3QoaXRlcnRvb2xzLmNvbWJpbmF0aW9ucyhyYW5nZShuKSwgcikpCiAgICAgIGFzc2VydCB3YW50ID09IGdvdCwgIm49e30gcj17fToge30gIT0ge30iLmZvcm1hdChuLCByLCB3YW50LCBnb3QpCgpzaG93X2FsbCg1LCAzKQp0ZXN0X2ZpcnN0KCkKdGVzdF9jb21iaW5hdGlvbigpCg==