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}$$