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

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.