trying https://github.com/keeganryan/flatter for faster lattice reduction than LLL

flatter.sage

from subprocess import check_output, DEVNULL, CalledProcessError
import itertools
import IPython
import time


def to_fplll_format(M):
    m, n = M.dimensions()
    ret = ""
    s = "["
    for i in range(m):
        s += "["
        for j in range(n):
            s += str(M[i, j])
            if j < n - 1:
                s += " "
        s += "]"
        ret += s + "\n"
        s = ""
    ret += "]"
    return ret


def from_fplll_format(s):
    rows = []
    for line in s.splitlines():
        line = line.lstrip("[").rstrip("\n").rstrip("]")
        if len(line) == 0:
            break

        row = [int(x) for x in line.split(" ") if len(x) > 0 and x != "]"]
        rows += [row]
    m = len(rows)
    n = len(rows[0])
    for row in rows:
        assert len(row) == n

    L = Matrix(ZZ, m, n)
    for i in range(m):
        for j in range(n):
            L[i, j] = rows[i][j]
    return L


def sort_by_norm(M):
    # idk why .norm() is really slow
    M2 = [(sum([y ^ 2 for y in x]), x) for x in M]
    M2.sort()
    return matrix([y for x, y in M2])


def LLL(M):
    return M.LLL()


def flatter(M):
    # compile https://github.com/keeganryan/flatter and put it in $PATH
    s = to_fplll_format(M)
    try:
        ret = check_output(["flatter"], input=s.encode(), stderr=DEVNULL)
        return from_fplll_format(ret.decode())
    except CalledProcessError:
        print("Flatter error, matrix written to /tmp/flatter_error")
        with open("/tmp/flatter_error", "w") as f:
            f.write(s)
        return M.change_ring(ZZ)


def subset_sum():
    print("=" * 40)
    print("Subset sum")
    print("=" * 40)
    # (very) low density subset sum
    n = 128
    B = vector([ZZ(randrange(0, 1 << 2048)) for _ in range(n)])
    s = random_vector(Zmod(2), n).change_ring(ZZ)
    C = B * s
    density = n / log(max(B), 2)
    print("density =", density.n())

    M = Matrix(ZZ, n + 1, n + 1)
    for i in range(n):
        M[i, i] = 2
        M[i, n] = B[i]
        M[n, i] = -1
    M[n, n] = -C
    M[:, n] *= M.det()

    print("flatter")
    print(IPython.get_ipython().run_line_magic("time", "flatter(M)[0]"))
    print("LLL")
    print(IPython.get_ipython().run_line_magic("time", "M.LLL()[0]"))
    print(vector([2 * x - 1 for x in s]))


def small_roots(f, bounds, m=1, d=None, reduction=LLL):
    # modified from https://github.com/defund/coppersmith/blob/master/coppersmith.sage
    if not d:
        d = f.degree()

    R = f.base_ring()
    N = R.cardinality()

    f /= f.coefficients().pop(0)
    f = f.change_ring(ZZ)

    G = Sequence([], f.parent())
    for i in range(m + 1):
        base = N ^ (m - i) * f ^ i
        for shifts in itertools.product(range(d), repeat=f.nvariables()):
            g = base * prod(map(power, f.variables(), shifts))
            G.append(g)

    B, monomials = G.coefficient_matrix()
    monomials = vector(monomials)

    factors = [monomial(*bounds) for monomial in monomials]
    for i, factor in enumerate(factors):
        B.rescale_col(i, factor)

    B = reduction(B.dense_matrix())

    B = B.change_ring(QQ)
    for i, factor in enumerate(factors):
        B.rescale_col(i, 1 / factor)

    H = Sequence([], f.parent().change_ring(QQ))
    for h in filter(None, B * monomials):
        H.append(h)
        I = H.ideal()
        if I.dimension() == -1:
            H.pop()
        elif I.dimension() == 0:
            roots = []
            for root in I.variety(ring=QQ, algorithm="msolve", proof=False):
                root = tuple(R(root[var]) for var in f.variables())
                roots.append(root)
            return roots

    return []


def univariate():
    print("=" * 40)
    print("Univariate Coppersmith")
    print("=" * 40)
    p = random_prime(2 ^ 1024, proof=False)
    q = random_prime(2 ^ 1024, proof=False)
    N = p * q
    bounds = (floor(N ^ 0.315),)
    roots = tuple(randrange(bound) for bound in bounds)
    R = Integers(N)
    P = PolynomialRing(R, 1, "x")
    x = P.gen()
    monomials = [x, x ^ 2, x ^ 3]
    f = sum(randrange(N) * monomial for monomial in monomials)
    f -= f(*roots)
    print("flatter")
    print(
        IPython.get_ipython().run_line_magic(
            "time", "small_roots(f, bounds, m=12, reduction=flatter)"
        )
    )
    print("LLL")
    print(IPython.get_ipython().run_line_magic("time", "small_roots(f, bounds, m=12)"))


def bivariate():
    def flatter_echelon(M):
        st = time.time()
        M = M.echelon_form(algorithm="pari", include_zero_rows=False)
        print("Echelon form done", time.time() - st)
        return flatter(M)

    def flatter_echelon_sort_by_norm(M):
        st = time.time()
        M = M.echelon_form(algorithm="pari", include_zero_rows=False)
        M = sort_by_norm(M)
        print("Echelon form + sort done", time.time() - st)
        return flatter(M)

    print("=" * 40)
    print("Bivariate Coppersmith")
    print("=" * 40)
    p = random_prime(2 ^ 1024, proof=False)
    q = random_prime(2 ^ 1024, proof=False)
    N = p * q
    bounds = (floor(N ^ 0.12), floor(N ^ 0.12))
    roots = tuple(randrange(bound) for bound in bounds)
    R = Integers(N)
    P = PolynomialRing(R, "x, y")
    x, y = P.gens()
    monomials = [x, y, x * y, x ^ 2, y ^ 2]
    f = sum(randrange(N) * monomial for monomial in monomials)
    f -= f(*roots)
    print("flatter")
    print(
        IPython.get_ipython().run_line_magic(
            "time", "small_roots(f, bounds, m=4, d=3, reduction=flatter)"
        )
    )
    print("flatter_echelon")
    print(
        IPython.get_ipython().run_line_magic(
            "time", "small_roots(f, bounds, m=4, d=3, reduction=flatter_echelon)"
        )
    )
    print("flatter_echelon_sort_by_norm")
    print(
        IPython.get_ipython().run_line_magic(
            "time",
            "small_roots(f, bounds, m=4, d=3, reduction=flatter_echelon_sort_by_norm)",
        )
    )
    print("LLL")
    print(
        IPython.get_ipython().run_line_magic("time", "small_roots(f, bounds, m=4, d=3)")
    )


if __name__ == "__main__":
    # copy this into sage repl or IPython wouldn't be available
    subset_sum()
    univariate()
    bivariate()

"""
========================================
Subset sum
========================================
density = 0.0625004965228275
flatter
CPU times: user 27.8 ms, sys: 196 µs, total: 28 ms
Wall time: 9.73 s
(1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, 1, -1, 1, 1, -1, -1, 1, -1, -1, -1, -1, -1, 1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, -1, -1, 1, -1, -1, 1, -1, 1, 1, -1, -1, -1, -1, 1, 1, -1, -1, -1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, 1, -1, 1, 1, -1, 1, -1, 1, -1, 1, 1, 1, 1, -1, 1, -1, -1, 1, 1, -1, 1, 1, 1, -1, 1, -1, -1, 1, -1, 1, 1, 1, -1, 1, -1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, -1, -1, 1, 1, 1, -1, 1, -1, -1, -1, 1, -1, -1, -1, -1, -1, 0)
LLL
CPU times: user 29.1 s, sys: 9.27 ms, total: 29.2 s
Wall time: 29.2 s
(-1, -1, 1, -1, 1, 1, -1, -1, 1, 1, 1, -1, 1, -1, -1, 1, 1, -1, 1, 1, 1, 1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, -1, -1, 1, 1, -1, 1, 1, 1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, -1, -1, -1, -1, 1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, -1, 1, -1, -1, -1, 1, -1, 1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, -1, 1, -1, 1, 1, 1, -1, 1, 1, 1, 1, 1, 0)
(-1, -1, 1, -1, 1, 1, -1, -1, 1, 1, 1, -1, 1, -1, -1, 1, 1, -1, 1, 1, 1, 1, 1, -1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, 1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, 1, -1, -1, 1, 1, -1, 1, 1, 1, -1, -1, -1, 1, -1, -1, 1, -1, 1, -1, 1, -1, -1, -1, -1, 1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, -1, 1, 1, -1, 1, -1, -1, -1, 1, -1, 1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, 1, -1, -1, -1, 1, -1, 1, 1, 1, -1, 1, 1, 1, 1, 1)
========================================
Univariate Coppersmith
========================================
flatter
CPU times: user 643 ms, sys: 20.1 ms, total: 663 ms
Wall time: 2.61 s
[(7135701385285258662705352079385971074031295450357866385772866275886697672493466130132019348479630137979490964611723360502755124684773373981646413594339672484490186402506928096443229421857274112,)]
LLL
CPU times: user 22.3 s, sys: 0 ns, total: 22.3 s
Wall time: 22.7 s
[(7135701385285258662705352079385971074031295450357866385772866275886697672493466130132019348479630137979490964611723360502755124684773373981646413594339672484490186402506928096443229421857274112,)]
========================================
Bivariate Coppersmith
========================================
flatter
Flatter error, matrix written to /tmp/flatter_error
CPU times: user 63.9 ms, sys: 0 ns, total: 63.9 ms
Wall time: 877 ms
[(0, 0)]
flatter_echelon
Echelon form done 0.053955078125
CPU times: user 186 ms, sys: 0 ns, total: 186 ms
Wall time: 20.9 s
[(16911423015429218599997249455414926104388466153079447100588523133119590144, 10859796556787345497472111002372524778473410645420854695590041467319954357)]
flatter_echelon_sort_by_norm
Echelon form + sort done 0.0610659122467041
CPU times: user 182 ms, sys: 9.98 ms, total: 192 ms
Wall time: 6.47 s
[(16911423015429218599997249455414926104388466153079447100588523133119590144, 10859796556787345497472111002372524778473410645420854695590041467319954357)]
LLL
CPU times: user 1.06 s, sys: 7 µs, total: 1.06 s
Wall time: 1.14 s
[(16911423015429218599997249455414926104388466153079447100588523133119590144, 10859796556787345497472111002372524778473410645420854695590041467319954357)]
"""