Translating a GAP-output for Sage

edit

This answer is threefold:

• first follow the instructions and stick to using GAP
• then propose a generator written in Sage
• finally compare the two generators

Call GAP from Sage via libGAP

We run the requested GAP commands via libGAP.

Tiny fixes and simplifications:

• degree $n$ polynomials have $n + 1$ coefficients
• remove some unneeded bits

Version (a) uses the code in the question with tiny simplifications:

def poly_list_a(n):
    r"""
    Return a list of ZZ-irreducible monic polynomials of this degree
    with coefficients 0, 1 or -1 and nonzero constant coefficient.
    """
    commands = f"""
    n := {n + 1}
    d := n - 1
    U := Tuples([-1, 1, 0], n)
    UU := Filtered(U, x -> (x[1] = 0) = false and x[n] = 1)
    W := []
    for i in UU do Append(W, [UnivariatePolynomial(Rationals, i)]); od
    WW := Filtered(W, x -> IsIrreducible(x) = true)
    """[1:-1]
    for command in commands.split('\n'):
        if command:
            _ = libgap.eval(command)
    return libgap.eval("WW").sage()

Version (b) tries to be a bit faster by not generating all tuples:

def poly_list_b(n):
    r"""
    Return a list of ZZ-irreducible monic polynomials of this degree
    with coefficients 0, 1 or -1 and nonzero constant coefficient.
    """
    commands = f"""
    n := {n}
    d := n - 1
    T := Tuples([-1, 1, 0], d)
    P := []
    for t in T do Append(P, Filtered(List([-1, 1], c -> UnivariatePolynomial(Integers, Concatenation([c], t, [1]))), IsIrreducible)); od
    """[1:-1]
    for command in commands.split('\n'):
        if command:
            _ = libgap.eval(command)
    return libgap.eval("P").sage()

Version (c) condenses the code quite a bit more:

def poly_list_c(n):
    r"""
    Return a list of ZZ-irreducible monic polynomials of this degree
    with coefficients 0, 1 or -1 and nonzero constant coefficient.
    """
    commands = f"""
    n := {n}
    T := Cartesian([[-1], [1]], Tuples([-1, 1, 0], n - 1), [[1]])
    P := Filtered(List(T, t -> UnivariatePolynomial(Integers, Concatenation(t))), IsIrreducible)
    """[1:-1]
    for command in commands.split('\n'):
        if command:
            _ = libgap.eval(command)
    return libgap.eval("P").sage()

Usage:

sage: poly_list = poly_list_c  # pick a or b or c

sage: print('\n'.join(str(p) for p in poly_list(1)))
x - 1
x + 1

sage: print('\n'.join(str(p) for p in poly_list(2)))
x^2 - x - 1
x^2 + x - 1
x^2 - x + 1
x^2 + 1
x^2 + x + 1

Sage code for a similar generator

Here we propose a generator written directly in Sage.

This calls other software too, but avoids mixing code in several languages at the user level.

We use the built-in generator for polynomials of given maximal degree over the finite field with three elements, and convert its polynomials to polynomials over ZZ with a little trick.

Transforming a polynomial to a dictionary skips 0 coefficients; the 1 and 2 coefficients can be mapped to 1 and -1 by the map x -> 3 - 2*x (after converting from GF(3) to ZZ).

def polys(n):
    r"""
    Generator for ZZ-irreducible monic polynomials of this degree
    with coefficients 0, 1 or -1 and nonzero constant coefficient.

    INPUT:

    - ``nn`` -- a positive integer (the degree)

    EXAMPLES::

        sage: print('\n'.join(str(p) for p in polys(1)))
        x + 1
        x - 1
        sage: print('\n'.join(str(p) for p in polys(2)))
        x^2 + 1
        x^2 - 1
        x^2 + x + 1
        x^2 + x - 1
        x^2 - x + 1
        x^2 - x - 1

    TESTS::

        sage: print('\n'.join(str(p) for p in polys(0)))
        Traceback (most recent call last):
        ...
        ValueError: Expected positive integer, got: 0
    """
    import numbers
    if not isinstance(n, numbers.Integral) or n < 1:
        raise ValueError(f"Expected positive integer, got: {n}")
    Zx = ZZ['x']
    for p in GF(3)['x'].polynomials(max_degree=n - 2):
        coeffs = {k + 1: 3 - 2*ZZ(c) for k, c in p.dict().items()}
        coeffs[n] = 1
        q = Zx(coeffs)
        for p in [q + 1, q - 1]:
            if p.is_irreducible():
                yield p

Usage:

sage: print('\n'.join(str(p) for p in polys(1)))
x + 1
x - 1

sage: print('\n'.join(str(p) for p in polys(2)))
x^2 + 1
x^2 + x + 1
x^2 + x - 1
x^2 - x + 1
x^2 - x - 1

Comparison

The good thing is both methods give the same polynomials:

sage: all(set(polys(n)) == set(poly_list(n)) for n in range(1, 8))
True

The poly_list_a function is a bit wasteful by removing tuples starting by 0 or not ending by 1 only after generating them, so it uses 2 out of 9 of the tuples it generates.

This is fixed in poly_list_b and poly_list_c. The gain is not spectacular since irreducibility tests dominate.

The serious drawback in poly_list_* is generating the whole list of polynomials rather than generating the polynomials one by one.

To test the speed of the various functions, we run through all the polynomials they give and do something fast with each polyomial, i.e. check that it is nonzero.

For n = 5, here are the timings:

sage: timeit('all(poly_list_a(5))')
5 loops, best of 3: 449 ms per loop

sage: timeit('all(poly_list_b(5))')
5 loops, best of 3: 420 ms per loop

sage: timeit('all(poly_list_c(5))')
5 loops, best of 3: 409 ms per loop

sage: timeit('all(polys(5))')
5 loops, best of 3: 55.4 ms per loop

The generator programmed in Sage is much faster; not sure what dominates in the speed difference:

• the overhead of calling libgap
• the overhead of storing a list
• other factors?

The GAP code might possibly be improved, e.g. by using GAP iterators but my GAP skills are limited.