Native support for solving linear systems over Laurent polynomial rings

[nzy1997]

Feature request

I would like to request native support in Oscar.jl for solving linear systems over Laurent polynomial rings, for example $R = k[x^{\pm 1}, y^{\pm 1}, z^{\pm 1}]$.
Given matrices $(A, b)$ with entries in such a ring, I am interested in solving $A \cdot u = b$ and, more generally, in computing the module of solutions over $R$.
Oscar has excellent support for linear algebra and module computations over polynomial rings, but it is currently unclear whether these operations are intended to work directly over Laurent polynomial rings.

I previously asked about this on Julia Discourse but did not find an answer: https://discourse.julialang.org/t/solving-a-u-b-over-a-laurent-polynomial-ring/134428

Current workaround

At the moment, I am able to solve such systems only via a fairly involved workaround, summarized below.
The idea is to:

1. Represent the Laurent polynomial ring $R$ as a quotient $P / I$ of a polynomial ring.
2. Map the linear system $A u = b$ to a module membership problem over $P$.
3. Enlarge the submodule generated by the columns of $A$ by adding generators corresponding to $I \cdot e_i$ for each basis vector.
4. Use coordinates to recover a solution and lift it back to $R$.
   A simplified version of the implementation is:

using Oscar

function solve_laurent_linear_equation(A, b)
    R = base_ring(A)
    m, n = size(A)

    f = Oscar._polyringquo(R)
    RQ = codomain(f)
    P = base_ring(RQ)
    I = modulus(RQ)
    
    A_quo = map_entries(f, A)
    A_poly = map_entries(x -> P(Oscar.lift(x)), A_quo)
    
    F = free_module(P, m)
    gens = [F(collect(A_poly[:, j])) for j in 1:n]
    
    for g in gens(I)
        for i in 1:m
            v = zeros(P, m)
            v[i] = g
            push!(gens, F(v))
        end
    end
    
    M, _ = sub(F, gens)
    
    b_quo = map_entries(f, b)
    b_poly = map_entries(x -> P(Oscar.lift(x)), b_quo)
    b_vec = F(collect(b_poly[:, 1]))
    
    coeffs = Oscar.coordinates(b_vec, M)
    u_poly = [coeffs[i] for i in 1:n]
    u_quo = [RQ(u_poly[i]) for i in 1:n]
    
    return matrix(R, n, 1, [preimage(f, u_quo[i]) for i in 1:n])
end

This works on small examples, e.g.:

F = GF(2)
R, (x, y, z) = laurent_polynomial_ring(F, ["x", "y", "z"])

A = matrix(R, 3, 3, [
    x       y^-1   1;
    z      x^-1   y;
    1      z^-1   x*y^-1
])

u = matrix(R, 3, 1, [x*y; z^-1; 1])
b = A * u

solve_laurent_linear_equation(A, b) == u

Issues with the workaround

While this approach demonstrates that the problem is mathematically manageable, it has several drawbacks:

• It relies on internal and undocumented APIs (_polyringquo, lift, preimage).
• It requires manually encoding quotient relations into a module, which is error-prone and non-intuitive.
• In experiments with larger matrices, the main performance bottleneck appears to be the repeated type conversions between rings, while the coordinates computation itself takes relatively little time.
  Overall, this feels like functionality that should be provided at the library level rather than implemented manually by users.

Request

It would be very helpful to have:

• native support for solving linear systems (or computing solution modules) over Laurent polynomial rings, or
• an official abstraction / recommended approach for handling such problems in Oscar.
  Thank you very much for considering this request, and for your work on Oscar.

[thofma]

I agree that we should have this functionality. I think an alternative (but equivalent) approach would be to add the solving for quotients of polynomial rings and then just use this in the Laurent case?