Function gfp_solve 

pub fn gfp_solve(
    f: NsField,
    equations: &[(Vec<(usize, u64)>, u64)],
    nvars: usize,
) -> Option<Vec<u64>>

Solve a linear system Σ aᵢ·x_{vᵢ} = rhs over f by incremental-echelon Gaussian elimination with modular-inverse pivot normalization. Each equation is a sparse coefficient vector over nvars unknowns plus its right-hand side. Returns a solution (free variables 0) or None on inconsistency. The field-general sibling of the GF(2) engine’s bit-packed gf2_solve: rows join an echelon basis one at a time (memory stays O(rank · nvars) however many equations stream through), each stored row is pivot-normalized with every entry strictly below its pivot column, so back-substitution in increasing pivot order reads a solution off directly.