Howgrave-Graham-Joux algorithm in Python
Subset sum in roughly O(2^(0.337n)) time for dense instances, using the representation technique.
- File
howgrave_graham_joux_algorithm.py, 76 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/howgrave-graham-joux-algorithm/. The authors are listed in prov/ZVD66Y.jsonl.
import itertools
import math
import random
def _list(
a: list[int], left: list[int], i: int, right: list[int], j: int, r: int, m: int
) -> list[int]:
"""Return the masks x with i ones in left, j ones in right, and sum(a[x]) = r (mod m)."""
g: dict[int, list[int]] = {} # @line L12
for z in itertools.combinations(right, j):
g.setdefault(sum(a[p] for p in z) % m, []).append(sum(1 << p for p in z))
out: list[int] = [] # @line L13
for x in itertools.combinations(left, i): # @line L14
mx = sum(1 << p for p in x)
for mz in g.get((r - sum(a[p] for p in x)) % m, []):
out.append(mx | mz)
return out # @line L15
def _weight_sum(a: list[int], mask: int) -> int:
return sum(a[p] for p in range(len(a)) if mask >> p & 1)
def howgrave_graham_joux(
a: list[int],
t: int,
w: int,
rounds: int | None = None,
rng: random.Random | None = None,
fallback: bool = True,
) -> list[int] | None:
"""Return the sorted positions of a w-subset of a with sum t, or None if there is none.
One level of representations: a solution e of weight w with
floor(w/2) ones in the left half is written as x + y, where x takes
floor(p/2) of its left ones and floor(q/2) of its right ones. A random
modular constraint on the sum of x keeps about one of these
representations. The answer is always correct: after `rounds` random
rounds, the search falls back to trying all w-subsets. The cost of a
round is the size of the lists it builds; the bounds on the page are
heuristic. With fallback=False the search is Monte Carlo: it may
return None although a solution exists, but a returned set is correct.
"""
n = len(a)
if w < 0 or w > n:
return None
rng = rng if rng is not None else random.Random()
rounds = 4 * (n + 1) if rounds is None else rounds
p, q = w // 2, w - w // 2 # @line L1
p1, q1 = p // 2, q // 2 # @line L2
m = math.comb(p, p1) * math.comb(q, q1) # @line L3
for _ in range(rounds): # @line L4
perm = list(range(n)) # @line L5
rng.shuffle(perm)
left, right = perm[: n // 2], perm[n // 2 :]
r = rng.randrange(m) # @line L6
lx = _list(a, left, p1, right, q1, r, m) # @line L7
ly = _list(a, left, p - p1, right, q - q1, (t - r) % m, m) # @line L8
h: dict[int, list[int]] = {} # @line L9
for y in ly:
h.setdefault(_weight_sum(a, y), []).append(y)
for x in lx: # @line L10
for y in h.get(t - _weight_sum(a, x), []):
if x & y == 0:
e = x | y
return [i for i in range(n) if e >> i & 1]
if not fallback:
return None
for c in itertools.combinations(range(n), w): # @line L11
if sum(a[i] for i in c) == t:
return list(c)
return None
test_howgrave_graham_joux_algorithm.py 105 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/howgrave-graham-joux-algorithm/. The authors are listed in prov/ZVD66Y.jsonl.
import itertools
import math
import random
from howgrave_graham_joux_algorithm import _list, howgrave_graham_joux
def brute_force(a: list[int], t: int, w: int) -> bool:
return any(sum(a[i] for i in c) == t for c in itertools.combinations(range(len(a)), w))
def is_valid(a: list[int], t: int, w: int, s: list[int] | None) -> bool:
return (
s is not None
and len(s) == w
and s == sorted(set(s))
and all(0 <= i < len(a) for i in s)
and sum(a[i] for i in s) == t
)
def test_example() -> None:
a = [31, 17, 45, 8, 26, 53, 12, 39]
for seed in range(20):
assert howgrave_graham_joux(a, 114, 4, rng=random.Random(seed)) == [0, 2, 4, 6]
assert howgrave_graham_joux(a, 102, 4, rng=random.Random(0)) is None
def test_example_lists() -> None:
a = [31, 17, 45, 8, 26, 53, 12, 39]
left, right = [0, 1, 2, 3], [4, 5, 6, 7]
sizes = [len(_list(a, left, 1, right, 1, r, 4)) for r in range(4)]
assert sizes == [4, 4, 4, 4]
sums = sorted(sum(a[p] for p in range(8) if x >> p & 1) for x in _list(a, left, 1, right, 1, 1, 4))
assert sums == [29, 57, 57, 61]
def test_edges() -> None:
assert howgrave_graham_joux([], 0, 0) == []
assert howgrave_graham_joux([], 0, 1) is None
assert howgrave_graham_joux([4], 4, 1) == [0]
assert howgrave_graham_joux([4], 0, 0) == []
assert howgrave_graham_joux([4], 4, 0) is None
assert howgrave_graham_joux([1, 2], 3, 3) is None
assert howgrave_graham_joux([1, 2], 3, -1) is None
same = [7] * 6
assert is_valid(same, 21, 3, howgrave_graham_joux(same, 21, 3, rng=random.Random(3)))
assert howgrave_graham_joux(same, 20, 3, rng=random.Random(3)) is None
big = [2**62, 2**62 - 1, 3, -(2**62), 2**61, 1, 5, 9]
t = 2**62 + 2**61 + 3 + 9
assert is_valid(big, t, 4, howgrave_graham_joux(big, t, 4, rng=random.Random(4)))
def test_list_against_brute_force() -> None:
rng = random.Random(77)
for _ in range(300):
n = rng.randrange(0, 9)
a = [rng.randrange(-20, 40) for _ in range(n)]
perm = list(range(n))
rng.shuffle(perm)
left, right = perm[: n // 2], perm[n // 2 :]
i, j = rng.randrange(0, 3), rng.randrange(0, 3)
m = rng.randrange(1, 7)
r = rng.randrange(m)
want = sorted(
sum(1 << p for p in x) | sum(1 << p for p in z)
for x in itertools.combinations(left, i)
for z in itertools.combinations(right, j)
if (sum(a[p] for p in x) + sum(a[p] for p in z) - r) % m == 0
)
assert sorted(_list(a, left, i, right, j, r, m)) == want
def test_against_brute_force() -> None:
rng = random.Random(20100530)
for _ in range(800):
n = rng.randrange(0, 11)
a = [rng.randrange(-30, 60) for _ in range(n)]
w = rng.randrange(0, n + 2)
t = rng.randrange(-40, 200)
got = howgrave_graham_joux(a, t, w, rng=rng)
if brute_force(a, t, w):
assert is_valid(a, t, w, got)
else:
assert got is None
def test_rounds_alone_find_planted_solutions() -> None:
rng = random.Random(5)
found = 0
trials = 200
for _ in range(trials):
n = 12
a = [rng.randrange(1, 2**20) for _ in range(n)]
e = rng.sample(range(n), n // 2)
t = sum(a[i] for i in e)
got = howgrave_graham_joux(a, t, n // 2, rounds=40, rng=rng, fallback=False)
if got is not None:
assert is_valid(a, t, n // 2, got)
found += 1
assert found >= trials * 9 // 10
assert math.comb(3, 1) * math.comb(3, 1) == 9
What the next level needs
- the tests did not run in this build