Solving A·u = b over a Laurent polynomial ring

General Usage
1

Hi, I am trying to solve linear equations over a Laurent polynomial ring in Oscar/Julia.

Given a matrix A and a right-hand side b over R = GF(p)[x₁^{±1}, …, xₙ^{±1}],
I want to compute a vector u such that A * u = b.

My current approach is:

  1. Convert the Laurent ring into a quotient ring R_Q = P/I using Oscar._polyringquo.
  2. Lift all entries to the polynomial ring P.
  3. Build a submodule M ⊆ P^m generated by the columns of A together with I * e_i.
  4. Use Oscar.coordinates(b, M) to express b in the Gröbner basis of M.
  5. Convert the solution coefficients back through
    P → R_Q → R.

This method works correctly (it reproduces the true Laurent solution), but the final step:

preimage(f, u_coeffs_quo[I])

is extremely slow and clearly the main bottleneck.

My questions:

  • Is there a more efficient way to convert quotient-ring elements back to Laurent polynomials?
  • Can I avoid calling preimage per coefficient, or use a faster map?
  • Is there a more direct method for solving A*x = b in a Laurent polynomial ring using Oscar/Singular, without going through heavy quotient-ring conversions?

Any advice or pointers to recommended techniques would be greatly appreciated!

Code

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_module = free_module(P, m)
    generators = [F_module(collect(A_poly[:, j])) for j in 1:n]
    
    ideal_gens = gens(I)
    for g in ideal_gens
        for i in 1:m
            v = zeros(P, m)
            v[i] = g
            push!(generators, F_module(v))
        end
    end
    
    M, inc = sub(F_module, generators)
    
    b_quo = map_entries(f, b)
    b_poly = map_entries(x -> P(Oscar.lift(x)), b_quo)
    b_vec = F_module(collect(b_poly[:, 1]))
    
    coeffs_vec = Oscar.coordinates(b_vec, M)
    u_coeffs_poly = [coeffs_vec[i] for i in 1:n]
    
    u_coeffs_quo = [RQ(u_coeffs_poly[i]) for i in 1:n]
    u_laurent = matrix(R, n, 1, [preimage(f, u_coeffs_quo[i]) for i in 1:n])

    return u_laurent
end

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

us = solve_laurent_linear_equation(A, b)
us == u