Laderman's algorithm in Python
Multiplies two 3×3 matrices with 23 scalar multiplications instead of 27; used recursively on blocks it costs O(n^(log3 23)), about n^2.854, which is worse than Strassen.
- File
laderman_algorithm.py, 131 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/laderman-algorithm/. The authors are listed in prov/JHP9MS.jsonl.
from typing import Optional
Matrix = list[list[int]]
# The 23 products M_r = (sum of c * A_ij) (sum of c * B_jk), with 1-based block indices.
# Each factor is a list of (coefficient, row, column).
PRODUCTS: list[tuple[list[tuple[int, int, int]], list[tuple[int, int, int]]]] = [
([(1, 1, 1), (1, 1, 2), (1, 1, 3), (-1, 2, 1), (-1, 2, 2), (-1, 3, 2), (-1, 3, 3)], [(1, 2, 2)]),
([(1, 1, 1), (-1, 2, 1)], [(-1, 1, 2), (1, 2, 2)]),
([(1, 2, 2)], [(-1, 1, 1), (1, 1, 2), (1, 2, 1), (-1, 2, 2), (-1, 2, 3), (-1, 3, 1), (1, 3, 3)]),
([(-1, 1, 1), (1, 2, 1), (1, 2, 2)], [(1, 1, 1), (-1, 1, 2), (1, 2, 2)]),
([(1, 2, 1), (1, 2, 2)], [(-1, 1, 1), (1, 1, 2)]),
([(1, 1, 1)], [(1, 1, 1)]),
([(-1, 1, 1), (1, 3, 1), (1, 3, 2)], [(1, 1, 1), (-1, 1, 3), (1, 2, 3)]),
([(-1, 1, 1), (1, 3, 1)], [(1, 1, 3), (-1, 2, 3)]),
([(1, 3, 1), (1, 3, 2)], [(-1, 1, 1), (1, 1, 3)]),
([(1, 1, 1), (1, 1, 2), (1, 1, 3), (-1, 2, 2), (-1, 2, 3), (-1, 3, 1), (-1, 3, 2)], [(1, 2, 3)]),
([(1, 3, 2)], [(-1, 1, 1), (1, 1, 3), (1, 2, 1), (-1, 2, 2), (-1, 2, 3), (-1, 3, 1), (1, 3, 2)]),
([(-1, 1, 3), (1, 3, 2), (1, 3, 3)], [(1, 2, 2), (1, 3, 1), (-1, 3, 2)]),
([(1, 1, 3), (-1, 3, 3)], [(1, 2, 2), (-1, 3, 2)]),
([(1, 1, 3)], [(1, 3, 1)]),
([(1, 3, 2), (1, 3, 3)], [(-1, 3, 1), (1, 3, 2)]),
([(-1, 1, 3), (1, 2, 2), (1, 2, 3)], [(1, 2, 3), (1, 3, 1), (-1, 3, 3)]),
([(1, 1, 3), (-1, 2, 3)], [(1, 2, 3), (-1, 3, 3)]),
([(1, 2, 2), (1, 2, 3)], [(-1, 3, 1), (1, 3, 3)]),
([(1, 1, 2)], [(1, 2, 1)]),
([(1, 2, 3)], [(1, 3, 2)]),
([(1, 2, 1)], [(1, 1, 3)]),
([(1, 3, 1)], [(1, 1, 2)]),
([(1, 3, 3)], [(1, 3, 3)]),
]
# C_ik is the sum of the products with these numbers (1-based); all coefficients are +1.
OUTPUTS: dict[tuple[int, int], list[int]] = {
(1, 1): [6, 14, 19],
(1, 2): [1, 4, 5, 6, 12, 14, 15],
(1, 3): [6, 7, 9, 10, 14, 16, 18],
(2, 1): [2, 3, 4, 6, 14, 16, 17],
(2, 2): [2, 4, 5, 6, 20],
(2, 3): [14, 16, 17, 18, 21],
(3, 1): [6, 7, 8, 11, 12, 13, 14],
(3, 2): [12, 13, 14, 15, 22],
(3, 3): [6, 7, 8, 9, 23],
}
def _classical(a: Matrix, b: Matrix, p: int, mod: Optional[int]) -> Matrix:
c = [[0] * p for _ in range(len(a))]
for i, row in enumerate(a):
ci = c[i]
for k, aik in enumerate(row):
if aik:
bk = b[k]
for j in range(p):
ci[j] += aik * bk[j]
if mod is not None:
c[i] = [x % mod for x in ci]
return c
def _combine(blocks: dict[tuple[int, int], Matrix], terms: list[tuple[int, int, int]],
rows: int, cols: int, mod: Optional[int]) -> Matrix:
out = [[0] * cols for _ in range(rows)]
for coef, i, j in terms:
blk = blocks[(i, j)]
for r in range(rows):
orow, brow = out[r], blk[r]
for s in range(cols):
orow[s] += coef * brow[s]
if mod is not None:
out = [[x % mod for x in row] for row in out]
return out
def laderman_algorithm(a: Matrix, b: Matrix, p: Optional[int] = None, mod: Optional[int] = None,
n0: int = 1) -> Matrix:
"""Return a * b by Laderman's 23-product recursion on 3 x 3 grids of blocks.
a is n x m and b is m x p (lists of rows); pass p when m = 0. The ring is
Z, or Z/mod Z when mod is set; commutativity is not used. Inputs with a
side of at most n0 (n0 >= 1) use the classical product; other sides are
padded with zeros to a multiple of 3 at each level. Cost: O(n^(log3 23))
ring operations for n x n inputs; recursion depth O(log max(n, m, p)).
"""
if n0 < 1:
raise ValueError("the cutoff n0 must be at least 1")
n, m = len(a), len(b)
if p is None:
p = len(b[0]) if m > 0 else 0
if any(len(row) != m for row in a) or any(len(row) != p for row in b):
raise ValueError("the shapes do not match")
return _laderman(a, b, n, m, p, mod, n0)
def _laderman(a: Matrix, b: Matrix, n: int, m: int, p: int, mod: Optional[int], n0: int) -> Matrix:
if min(n, m, p) <= n0: # @line L1
return _classical(a, b, p, mod)
if n % 3 or m % 3 or p % 3: # @line L2
n3, m3, p3 = -(-n // 3) * 3, -(-m // 3) * 3, -(-p // 3) * 3
ap = [row + [0] * (m3 - m) for row in a] + [[0] * m3 for _ in range(n3 - n)]
bp = [row + [0] * (p3 - p) for row in b] + [[0] * p3 for _ in range(m3 - m)]
cp = _laderman(ap, bp, n3, m3, p3, mod, n0)
return [row[:p] for row in cp[:n]]
h, k, q = n // 3, m // 3, p // 3 # @line L3
ab = {(i + 1, j + 1): [row[j * k:(j + 1) * k] for row in a[i * h:(i + 1) * h]]
for i in range(3) for j in range(3)}
bb = {(i + 1, j + 1): [row[j * q:(j + 1) * q] for row in b[i * k:(i + 1) * k]]
for i in range(3) for j in range(3)}
prods = []
for left, right in PRODUCTS: # @line L4
x = _combine(ab, left, h, k, mod)
y = _combine(bb, right, k, q, mod)
prods.append(_laderman(x, y, h, k, q, mod, n0))
cb = {}
for (i, j), idx in OUTPUTS.items(): # @line L5
cb[(i, j)] = _combine({(r, 0): prods[r - 1] for r in idx}, [(1, r, 0) for r in idx], h, q, mod)
out = [] # @line L6
for i in range(1, 4):
for r in range(h):
out.append(cb[(i, 1)][r] + cb[(i, 2)][r] + cb[(i, 3)][r])
return out
def count_multiplications(n: int, n0: int = 1) -> int:
"""Scalar multiplications for n x n inputs with n = 3^t * s, s <= n0: 23 T(n/3), base n^3."""
if n <= n0:
return n ** 3
return 23 * count_multiplications(n // 3, n0)
test_laderman_algorithm.py 87 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/laderman-algorithm/. The authors are listed in prov/JHP9MS.jsonl.
import random
import pytest
from laderman_algorithm import OUTPUTS, PRODUCTS, count_multiplications, laderman_algorithm as mul
def schoolbook(a, b, p, mod=None):
n, m = len(a), len(b)
c = [[0] * p for _ in range(n)]
for i in range(n):
for j in range(p):
s = 0
for k in range(m):
s += a[i][k] * b[k][j]
c[i][j] = s % mod if mod is not None else s
return c
def rand(rng, r, c, lo=-9, hi=9):
return [[rng.randint(lo, hi) for _ in range(c)] for _ in range(r)]
def test_example() -> None:
a = [[1, 2, 3], [4, 5, 6], [7, 8, 10]]
b = [[1, 0, 2], [0, 1, 1], [1, 1, 0]]
assert mul(a, b) == [[4, 5, 4], [10, 11, 13], [17, 18, 22]]
def test_scheme_shape() -> None:
assert len(PRODUCTS) == 23
adds = sum(len(u) - 1 + len(v) - 1 for u, v in PRODUCTS)
adds += sum(len(v) - 1 for v in OUTPUTS.values())
assert adds == 98
assert count_multiplications(3) == 23
assert count_multiplications(9) == 529
assert count_multiplications(27, n0=3) == 23 * 23 * 27
def test_edges() -> None:
assert mul([], []) == []
assert mul([[5]], [[-3]]) == [[-15]]
assert mul([[], []], [], p=2) == [[0, 0], [0, 0]]
with pytest.raises(ValueError):
mul([[1, 2]], [[1, 2]])
with pytest.raises(ValueError):
mul([[1]], [[1]], n0=0)
def test_identity_and_zero() -> None:
rng = random.Random(3)
for n in (2, 3, 4, 6, 9, 10):
a = rand(rng, n, n)
eye = [[int(i == j) for j in range(n)] for i in range(n)]
zero = [[0] * n for _ in range(n)]
assert mul(a, eye) == a
assert mul(eye, a) == a
assert mul(zero, a) == zero
def test_random_3x3_over_z() -> None:
rng = random.Random(77)
for _ in range(600):
a, b = rand(rng, 3, 3, -100, 100), rand(rng, 3, 3, -100, 100)
assert mul(a, b) == schoolbook(a, b, 3)
def test_random_shapes_and_recursion() -> None:
rng = random.Random(78)
for _ in range(80):
n, m, p = rng.randint(1, 11), rng.randint(1, 11), rng.randint(1, 11)
a, b = rand(rng, n, m), rand(rng, m, p)
assert mul(a, b, n0=rng.choice([1, 2, 3])) == schoolbook(a, b, p)
a, b = rand(rng, 27, 27), rand(rng, 27, 27)
assert mul(a, b, n0=3) == schoolbook(a, b, 27)
def test_mod_small_primes() -> None:
rng = random.Random(79)
for q in (2, 3, 5, 7):
for _ in range(40):
n = rng.choice([3, 4, 5, 6, 9])
a, b = rand(rng, n, n, 0, q - 1), rand(rng, n, n, 0, q - 1)
assert mul(a, b, mod=q) == schoolbook(a, b, n, q)
What the next level needs
- the tests did not run in this build