Solving a polynomial congruence with rational number unknowns for absolute factorisation I am implementing Gao's factorisation algorithm for bivariate rational polynomials $f\in\mathbb Q[x,y]$. An overview and the reference to the paper describing the algorithm are in this answer. I see value in the algorithm because it performs absolute factorisation – if the polynomial splits over some algebraic field, the algorithm will calculate it; I do not need to guess.
I am following the original paper closely and there is a step I am unable to implement explicitly (using SymPy).

Theorem 2.8. Suppose that $g_1,\dots,g_r$ form a basis for $G$ over $\mathbb F$ [which is $\mathbb Q$ in this question's context]. For
  any $g\in G$, there is a unique $r×r$ matrix $A=(a_{ij})$ over $\mathbb F$ such that
  $$gg_i\equiv\sum_{j=1}^ra_{ij}g_jf_x\mod f\tag1$$

$r$ is the number of absolutely irreducible factors of $f$. I have successfully implemented procedures to compute the $g_i$ (which arise as the nullspace of a linear system), and $g$ is a randomly chosen linear combination of the $g_i$. If $g$ is such that $A$'s characteristic polynomial $c_A(x)$ has no repeated roots, then it is shown that $f$ splits over $\mathbb Q(\alpha)$ where $c_A(\alpha)=0$.

What is the procedure to compute the $a_{ij}$ in $(1)$ when given $f$, the $g_i$ and the chosen $g$?

I believe the main difficulty is ensuring that the $a_{ij}$ are in $\mathbb Q$ – the routines I've examined in SymPy for Bézout decompositions of multivariate polynomials don't seem to be able to enforce this. The $\bmod f$ is also tripping me up.
There is a worked example given which may help with the explanation, with $f=9+23y^2+13yx^2+6y+7y^3+13y^2x^2+x^4+6yx^4+x^6$. This polynomial has three absolutely irreducible factors ($r=3$) with computed $g_i$
$$g_1=-12x-8xy-19xy^2-12x^3y-2x^5+x^3$$
$$g_2=12x+10xy+18xy^2+12x^3y+2x^5$$
$$g_3=-18x-12xy-22xy^2-14x^3y-2x^5$$
$$g=g_1+g_2=2xy-xy^2+x^3$$
The computed $A$ is
$$\begin{bmatrix}
-62/247&63/988&189/988\\
63/247&-17/247&-51/247\\
-54/247&135/494&79/247\end{bmatrix}$$
 A: The problem is actually fairly simple if the operations are performed in the right order. When $gg_i$ and the $g_jf_x$ polynomials are taken modulo $f$ first, the remainders' monomials will have degrees of the same order (rem(f,g) in SymPy), so that a linear system can be set up to find the $a_{ij}$. To illustrate, for the example's polynomials, the entries of $A$'s first row are the solutions to the linear system starting with
$$\begin{bmatrix}
0&12&-4\\
106&-108&196\\
96&-128&200\\
\vdots&\vdots&\vdots\end{bmatrix}\begin{bmatrix}a_{11}\\a_{12}\\a_{13}\end{bmatrix}=\begin{bmatrix}0\\4\\6\\\vdots\end{bmatrix}$$
where the columns from left to right are the reduced $g_jf_x$ and $gg_1$ polynomials respectively, and the displayed rows correspond to the $x^4y^3,x^4y^2,x^4y$ coefficients. Once that hurdle was cleared, I managed to complete the implementation and successfully find splitting fields for small bivariate polynomials.

Immediately afterwards, however, I saw an improvement from Jürgen Gerhard mentioned in the same paper that saves the hassle of finding a basis for the initial linear system $G$'s nullspace, guessing $g$ and constructing $A$ – the hassle that led me to ask this question in the first place. It involves taking any non-trivial $g$ in $G$ and computing the resultant $\operatorname{Res}_x(f,g-zf_x)$, from which the number of absolutely irreducible factors and the splitting field can be derived. The sheer size of the $G$ matrices I encountered with slightly larger polynomials also compelled me to write the implementation in PARI/GP instead, which managed to find a quintic splitting field for a degree-$25$ test polynomial I concocted. The PARI/GP implementation is available as gao.gp here.
